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1.  INTRODUCTION 


In  1985,  Mumford  and  Shah  introduced  a  mathematical  model  for  solving  the 
segmentation  problem  (Reference  1).  Segmentation  and  boundary  detection  algorithms 
are  basic  tools  for  extracting  global  features  out  of  digitized  data. 

The  Mumford  and  Shah  (M&Sh)  functional  (References  1  and  2)  has  the  form 


E(u,K)  =  a  J  ||u  -  gfdju  +  J||Vwf  dju  +  A  •  £(K) , 

n-K  n-K 


(i) 


where  Q  denotes  a  rectangle  in  Rd ,  d  e  {1,2,3, •••}  is  the  dimension  of  Q,  g:Q-»  RL 
denotes  the  image  (g  is  ju  -measurable),  c  e{l,2,3,---}  is  the  number  of  channels,  and 
u :  Q  Rc  is  a  smooth  approximation  to  g .  The  rectangle  Q  is  decomposed  into  a 
finite  collection  of  disjoint  open  sets  On  (n  =  1,  2,  N)  and  their  boundaries 


/v 

K  —  [J  dOn  ;  so  that  Q 


«=i 


U  o. 


.«=! 


u  K,  emd  K  is  sufficiently  nice  to  have  a  length. 


denoted  by  £(K) .  The  disjoint  open  sets  On  («=1,  2,  ...,  N)  and  their  boundaries 

N 

K  =  [J  dOn  together  with  u  are  called  a  segmentation  of  (  Q ,  g  ). 

M=1 


The  goal  is  to  construct  two  things: 

1 .  A  smoothed  ideal  image  u :  Q  -»  Rc 

2.  A  set  of  boundaries  K  c  Q 

u  and  K  are  found  by  minimizing  the  functional  E  . 
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The  first  term  on  the  right-hand  side  (RHS)  of  Equation  1  ensures  that  u  is  a  faithful 
representation  of  g ,  the  second  term  ensures  that  u  is  as  smooth  as  possible  on  each 
open  set  0„  (the  segments),  and  the  last  term  prevents  the  boundaries  from  growing  too 
large.  The  technique  is  multi-scale.  The  parameters  a  and  X  are  weighting  factors  that 
control  the  quality  of  the  approximation  and  the  coarseness  of  the  segmentation.  The 
technique  allows  the  extraction  of  features  at  different  levels  of  detail  (scale).  The 
technique  is  also  multichannel;  u  and  g  may  be  vector  valued  functions.  It  can  be  used 
to  segment  images  of  a  scene  when  registered  multiple  data  channels  for  the  same  scene 
are  available.  These  may  be  data  channels  from  various  sensors,  hue  channels,  or 
preprocessing  channels  such  as  wavelet  or  other  transform  channels.  Thus,  the  M&Sh 
model  captures  all  the  essential  features  that  must  be  considered. 

The  aim  of  this  report  is  to  present  an  implementation  of  the  region  growing  (RG) 
method  introduced  in  References  3  and  4  for  minimizing  the  simplified  M&Sh  functional 

Q-K 


Their  model  assumes  that  «:Q-»  Rc  is  a  piecewise  constant  (PC)  function  defined 

N 

on  (J  On  ,  constant  on  each  Oit ;  u\0„  ss  An ,  where  A„  is  the  average  value  of  g  on  On . 
/>=! 

Thus,  the  second  term  in  Equation  1  vanishes,  leading  to  Equation  2. 

The  idea  behind  the  RG  method  for  minimizing  Equation  2  is  to  select  two  adjacent 
regions  O, ,  Oj  and  merge  them;  that  is,  the  new  average  Ay  of  g  on  O,  u  0}  is  found; 
the  boundary  between  (9,  and  Oj}  denoted  by  d{0,,0i'),  is  eliminated;  and  the 
functional  E  is  evaluated.  If  E  decreases,  then  O,  and  Oj  are  permanently  merged, 
leading  to  a  coarser  segmentation.  Otherwise,  another  pair  is  considered.  Starting  with  the 
trivial  segmentation,  i.e.,  one  where  every  region  consists  of  only  one  point  (pixel),  this 
procedure  arrives  at  a  segmentation  with  minimal  (local)  energy  for  the  chosen  scale 
parameter  X .  If  X  is  large,  the  emphasis  is  on  the  length  of  the  boundaries  K  and  the 
minimization  of  E  leads  to  fewer  boundaries  and  hence  to  a  coarser  segmentation. 
Conversely,  a  small  X  puts  the  emphasis  on  the  first  term  of  Equation  2,  leading  to  an 
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approximation  u  that  follows  g  much  more  closely  and  allows  for  more  boundaries.  Thus, 
a  small  X  leads  to  a  finer  segmentation.  In  the  limit  when  X  is  zero,  the  minimization  of 
E  forces  u  =  g,  leading  to  the  finest  possible  segmentation,  the  trivial  segmentation  where 
every  region  consists  of  a  single  point  (pixel). 

The  merging  criterion  for  merging  two  adjacent  regions  Oj  and  is  defined  to  be 
the  difference  between  E(u,K)  and  the  new  value  E(u,K- d{0,,0 /))  after  a  merger. 

( Ay  on  Oj  u  O, 

If  u  =  An  on  each  On  and  u=\  '  ,  where  Au  is  the  average  value  of  g 

[w  otherwise 

on  Oj  u  Oj ,  then  the  merging  criterion  for  the  PC  case  (References  3  and  5)  is  given  by 

H(0j)n(0.)  ,,  .2 

M „  .  E(u, K)  -  £(». K  -  m, Oj ))  =  A •  t(S(0„ Oj ))  -  '  /’  ■  U,  -  At  \  ,  (3) 

where  ju(0,  )  denotes  the  area  of  the  open  set  O ) .  If  >0  there  is  a  decrease  in  E  and 

the  regions  Ot,  Ot  are  merged. 

While  the  expression  on  the  RHS  of  Equation  3  for  appears  in  Reference  3,  the 

authors  do  not  provide  a  derivation  of  it.  A  derivation  of  it  cannot  be  found  in  Reference 
4  either.  Reference  5  presents  a  complete  derivation  of  Equation  3  and  its  generalization 
to  the  case  where  the  approximations  un  are  not  PC;  that  is,  the  approximations  un  are 
restrictions  to  0„  of  functions  belonging  to  a  class  of  functions  defined  as  the  linear  span 
of  a  finite  set  of  judiciously  chosen  functions  /j,/2,  ...  ,f?  mapping  Q.  into  Rc .  For 
example,  this  can  include  PC,  piecewise  affine  (PA),  piecewise  polynomial,  piecewise 
exponential,  piecewise  sinusoidal  approximations. 

I 

In  the  general  setting  of  Reference  5,  each  un  =  u\On  has  the  form  un  =  yi  ank  fk  •  ^ 

*= 1 

we  let  An  =[«„,  an2  •••  ane]r  e Re  and  F  =  [/,:/2 :•••:/,],  then  un  can  be  written  as 
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un  =  FAn .  The  weighting  factor 
weighting  matrix  Hy ,  defined  by 


mMQj) 


appearing  in  Equation  3  becomes  a 


(4) 


where 


M„  =  J  Fr  Fd/j.  (n= 

On 


(5) 


The  merging  criterion  in  this  general  setting  has  the  form 

*E{u,K)-E(u,K-XOl,Oj))  =  A-iWOl,0J))-\Al-AJ(H'  ,  (6) 

where  the  weighted  norm  ||  -||w  is  defined  as  ||x||^  =xTHx  for  x  e  Re . 

Reference  5  also  includes  some  examples  showing  the  form  that  the  functions 
un  =  FAn  and  the  matrix  F  =  [/,:/2 :•••:/<]  take  in  various  settings.  These  include  the 

multichannel  PC,  the  multichannel  PA,  and  the  multichannel  piecewise  quadratic  (PQ) 
settings  in  various  dimensions. 

Here  we  use  the  results  of  Reference  5  and  illustrate  how  one  can  implement  the  RG 
method  in  dimension  2  for  PC  approximations  un  in  the  single-  and  multichannel 
settings.  The  multichannel  PA  implementation  is  described  in  Reference  6.  Sufficient 
information  can  be  found  in  Reference  6  to  implement  the  multichannel  PQ  algorithm 
following  similar  steps. 

In  Section  2  we  state  the  results  of  Reference  5  that  will  be  needed  here.  A  tiling 
problem  in  the  plane  that  arose  as  a  result  of  trying  to  speed  up  our  implementation  is 
discussed  in  Section  3. 
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A  general  description  of  the  program  structure  is  given  in  Section  4.  This  section  is 
generic  in  the  sense  that  all  of  the  different  versions  of  the  RG  method  that  we  have 
implemented  have  the  same  structure.  We  have  implemented  the  RG  method  in  the  PC 
single-  and  multichannel  settings  in  dimension  2  and  in  the  PA  single-  and  multichannel 
settings  in  dimensions  2  and  3. 

In  Section  5  we  describe  in  detail  the  Matlab  functions  that  implement  the  single¬ 
channel  PC-RG  method  in  two  dimensions  and  we  indicate  the  changes  involved  in  the 
multichannel  implementation.  These  changes  are  minimal.  Listings  of  the  functions  are 
included  in  the  Appendix  for  reference  and  for  completeness.  The  changes  involved  in  the 
multichannel  implementation  are  also  listed.  Thus,  this  report  serves  as  documentation  for 
the  two  dimensional  PC-RG  programs  in  the  single-  and  multichannel  settings. 


2.  PIECEWISE  APPROXIMATIONS 


In  this  section  we  show  how  to  obtain  the  piecewise  approximations  u  on  the  sets  On 
and  how  to  compute  the  approximation  Uy  on  a  union  0,  u  0}  from  the  approximations 
on  Oj  and  0} . 


For  each  n  =1 ,  2,...,  N,  let 

V,  =  \Frgdu,  (7) 

o„ 

where  ju  is  a  measure  on  Q  (e.g.,  Lebesgue  measure).  Fix  i  and  j  in  {1,  2,  N}  and 
define  u(j: Q  Rc  by  uij  ~  bkfk  .  Note  that  if  B  ~  [bx  b2  •••  b£]T ,  then  ul}  =  FB . 

k  =  ] 


N 

Define  u:  [J  On  — >  Rc  by  u  =  un  on  On  for  all  n. 

«=i 
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n  \  uv  on  O'  u  Oj 

Define  u:\  \0  ->  Rc  by  u  =  s 

,^1  [u  otherwise. 

Proposition  1.  (Piecewise  approximations  to  g ;  The  coefficient  vectors  An  and  B .) 

e 

(a)  If  u„  =  Y,ankfk  approximates  g  on  0„  in  the  sense  of  least  mean  square  error, 

1 

that  is,  if  {anl,aia,---,allf}  minimize  J||w„  -gf  dp,  then  the  vector  of  coefficients  An 

o„ 

satisfies  the  linear  system  M„An  =V„  (n  =  1,  2,...,  N ). 

(b)  Similarly,  if  B  minimizes  J  |wv  -gf  dp,  then  B  satisfies  the  linear  system 

OjVOj 

(Mi  +  MJ)B  =  (Vi+VJ). 

(c)  B  =  (Mi  +  Mj)'\MiAi  +  MjAj)  provided  the  inverse  exists. 

Proof.  (  See  Reference  5.) 

Proposition  2.  (The  merging  criterion.) 

If  u  and  u  are  as  above  and  E(u,K)=  J  ||w  -  gf  dp  +  X  •  t{K) ,  then  Equation  6 

n-K 

holds. 

Proof.  (  See  Reference  5.) 

2.1  MULTICHANNEL  PIECEWISE  CONSTANT  APPROXIMATIONS 

Suppose  u  is  a  PC  approximation  of  g;  that  is,  u  is  constant  on  each  set  0„  and  there 
are  c  channels.  In  this  case  F=  I  ,  the  cxc  identity  matrix;  because  if 
uH  =  [o„,  all2  ■  •  ■  anc]T ,  a  constant  vector  with  c  components,  then 
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T 

'0 

'O' 

0 

1 

0 

JZ 

II 

0 

+  a„2 

0 

+  •••  +anc 

0 

_0_ 

1 

—  an\f\  +  anlfl  "*  ancfc 


'0" 

'1 

0 

...  o' 

Therefore,  fk  = 

1 

<-  (k,h  -  position ) ,  F  =  [f1:f2:- 

0 

1 

...  0 

0 

_0 

0 

... 

,  and 


F‘  F  =  I . 


Since  Mn  =  J  Id/u  =  ju(0„)  •  /  for  all  n  ,  the  weight  matrix  Hy  is  a  multiple  of  the 


identity  matrix: 


Hu  = 


I  and  lUII2  = 


fKO,)+f<Oj) 


/KO,)+fKOj) 

The  merging  criterion  for  the  multichannel  PC  case  becomes 


*  • 


rtoj  +  AOj) 


Aj  A:  , 


where  A,  ( Aj  respectively)  is  the  vector  of  values  of  u  on  Oi  ( Oy  respectively).  Thus, 
Equation  3  is  a  special  case  of  Equation  6. 

By  Proposition  1,  A„  satisfies  ju(On )  ■  An  =  J  gdpi .  Thus,  A„  =  — - - \gdju  is  the 

o„  m,)  o„ 

average  value  of  g  on  On .  Similarly,  B  = - - - \gdju  is  the  average  value 

fKO,)+KOj)  0ii0j 
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of  g  on  O;  u  Oj .  In  this  case,  part  (c)  of  Proposition  1  simply  shows  how  to  obtain  the 
average  of  g  on  O,  u  Oj  from  the  average  ofg  on  0,  and  0} . 

As  an  example,  we  include  here  the  single  channel  PA  case  in  dimension  2  for 
comparison  and  to  emphasize  that  Equation  6  is  indeed  necessary  beyond  the  PC  setting 
and  that  Equation  3  applies  only  in  the  PC  setting. 


2.2.  SINGLE-CHANNEL  PIECEWISE  AFFINE  APPROXIMATIONS  (d  =  2) 


If  d  =  2,  then  O  is  a  rectangle  in  R2 .  Suppose  g:Q-^  R  and  u  is  a  PA 
approximation  to  g.  Then 


u„  (x,y)  =  a„,  •  1  +  a„2  ■  x  +  all3  •  .y  =  [1  *  y] 


a 


Ml 


a 


m2 


La«3j 


=  F{x,y)An . 


Therefore,  /,  (x,y)  =  1 ,  a  constant  function,  /2  (x,y)  =  x ,  and  /3(x,y)  =  y ,  (x,y)  e  Q . 
The  matrix  F(x,y)  =  [  1  x  y  ]  is  a  1  x  3  matrix,  and 


T 

1  x  y 

Fr(x,y)F(x,y)  = 

X 

_y_ 

[1  x  y]  = 

x  x2  xy 

y  yx  y\ 

(x ,y)  eQ. 


J  F  ‘\x,y)F{x,y)dxdy 
o„ 


x  y 
x2  xy 

yx  y2 


dxdy 
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1 1  dxdy  |  x  dxdy  J  y  dxdy 
o,  •>  0„  o„ 

=  |  x  dxdy  J  x1  dxdy  J  xy  dxdy  , 
o„  on  o„ 

|  y  dxdy  J  xy  dxdy  J  y2  dxdy 

0„  0„  O. 

so  Mn  is  no  longer  a  multiple  of  the  identity  matrix. 


K  =  f^T(x,y)g(x,y)dxdy  = 

On 


Jl  -g(x,y)dxdy 
On 

\x-g(x,y)dxdy 

o„ 

\y-gix,y)dxdy 

0, 


The  coefficient  vector  A„ 


a 


Hi 


a 


h2 


a 


n  3 


satisfies  the  linear  system  M„  A„  -  V„ . 


Since  the  matrices  M„  and  Hy  are  no  longer  multiples  of  the  identity  matrix, 

Equation  3  no  longer  applies  and  one  must  use  the  general  merging  criterion  in  Equation 
6.  (See  Reference  6  for  more  details  on  the  multichannel  PA  setting.) 


3.  TILING  PROBLEM 


When  implementing  the  RG  method,  the  question  of  which  pair  of  adjacent  regions 
Ot,Oj  should  be  considered  first  arises.  What  procedure  should  be  used  to  select  the 

sequence  of  pairs  of  adjacent  regions?  In  our  initial  implementation  we  decided  to  merge 
a  pair  of  regions  0,,0f  with  maximal  criterion  M  y .  This  procedure  leads  to  a  steepest- 
descent  type  of  procedure  because  the  maximal  criterion  My  represents  the  largest 
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decrease  in  E  when  two  regions  are  merged.  We  discovered,  however,  that  this  approach 
often  leads  to  regions  that  grow  too  large  and  with  too  many  small  neighboring  regions. 
When  regions  become  such  that  they  have  too  many  neighbors,  the  procedure  slows  down 
considerably.  We  concluded  that,  at  least  initially,  one  should  avoid  large  regions  with 
too  many  small  neighbors.  We  opted  to  use  a  seeding  procedure  at  the  initial  stages  of  the 
two-dimensional  merging  process.  The  seeding  procedure  (tiling)  is  described  next. 

The  RG  method  starts  with  an  initial  segmentation  of  Q .  In  the  discrete  PC  case,  the 
initial  segmentation  can  be  chosen  to  be  the  trivial  segmentation,  that  is,  one  in  which 
each  On  consists  of  one  single  pixel  (PC  case  only).  Then  each  pixel  has  four  neighbors: 
up,  down,  right,  and  left  neighbor.  Note  that  the  diagonal  neighbors  of  a  pixel  x  share  a 
boundary  with  x  (a  comer)  of  zero  length.  Consequently,  these  neighbors  will  never 
merge  with  x,  because  the  merging  criterion  in  this  case  is  never  positive  (see  Equation 
6). 


When  a  pixel  x  merges  with  its  four  neighbors,  it  forms  a  cross-like  region  (Figure  1) 
which  we  will  call  a  first  level  tile  (FL-tile).  The  pixel  x  is  the  center  of  the  tile. 


12 


NAWCWD  TP  8525 


3.1  ITERATED  TILING  PROBLEM  IN  TWO  DIMENSIONS 

The  tiling  problem  for  fici!2  can  be  stated  as  follows.  (Note  that  Q  may  be 
replaced  by  R2 .) 

1 .  Determine  the  centers  of  a  collection  of  FL-tiles  so  that  the  collection  of  FL-tiles 
cover  Q  and  do  not  overlap  (i.e.,  intersect  only  at  the  boundary).  Such  a  covering  will  be 
called  a  first-level  tiling  of  Q . 

2.  Determine  the  centers  of  a  collection  of  the  resulting  FL-tiles  so  that  the  second- 
level  tiles  (SL-tiles),  that  is,  a  FL-tile  merged  with  its  adjacent  FL-tile  neighbors,  cover 
Q  and  do  not  overlap.  Such  a  covering  will  be  called  a  second-level  tiling  of  Q . 

3.  Determine  how  to  iterate  the  tiling  procedure  implied  by  1  and  2  to  obtain  higher- 
level  tilings. 


3.2  ITERATED  TILING  PROBLEM  SOLUTION  IN  TWO  DIMENSIONS 

1 .  First-level  tiling  (discrete  case). 

Let  Q  =  { (ij)-pixels  :  0  <  i  <  r,  and  0  < j  <  s  }.  Let 
j3/(mod5)  for  i  =  0,1, 2, 3, 4 

a‘  Ws  f°r  1  ^ 5 

then,  the  centers  for  the  FL-tiles  can  be  chosen  to  be 


C,  =  { cjk  =  (i,ai  +  5k)  :i  s  Z ,k  e  Z  }  n  Q . 


Note  that  the  centers  in  C,  are  distributed  with  density  1/5  per  unit  area  and  the  area 
of  a  FL-tile  is  5.  Trivially  one  can  check  that  FL-tiles  with  centers  in  C,  cover  Q  (except 
for  edges)  and  do  not  overlap. 
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2.  Second-level  tiling. 

The  centers  C2  for  the  second  level  tiling  must  be  chosen  from  the  set  Cx .  There  are 
several  choices.  They  must  be  spaced  a  distance  five  in  both  directions.  For  example: 

C2={cik=(5i,5k):i  eZ,k  eZ}nCl. 

Remark  3.1.  The  process  iterates  to  higher-level  tiles.  The  centers  in  C2  form  a  square 
lattice  just  as  the  centers  of  the  initial  pixels  do.  The  second-level  tiles  are  square  like 
(Figure  2).  The  third  level  tiles  will  form  a  cross-like  region  (Figure  3).  Thus,  the  process 
clearly  iterates. 


FIGURE  2.  Second-Level  Tiles. 
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Remark  3.2.  The  goal  of  the  tiling  problem  is  to  prevent  regions  from  having  too  many 
small  neighbors.  In  the  tiling  process  previously  described,  the  regions  maintain  four,  and 
only  four,  neighbors  throughout  the  process.  Furthermore,  the  area  of  FL-tiles  is  5,  the 
SL-tiles  have  area  25,  and  the  k-level  tiles  have  area  5* .  Thus,  the  tile  area  grows 
exponentially  while  the  number  of  neighbors  remains  constant.  This  is  the  case  for  the 
worst-case  images,  which  ironically  for  the  RG  method  are  the  blank  images  or 
homogeneous  images  (images  with  constant  g)  that  force  the  method  to  merge  all  the 
initial  regions  into  one.  For  non-blank  images  the  iterated  tiling  process  will  be  corrupted 
when  lower-level  tiles  do  not  get  merged  into  higher-level  tiles  because  they  are  not 
supposed  to  (merging  criterion  <  0).  After  a  few  iterations  of  the  tiling  process  our 
algorithm  starts  merging  regions  with  maximal  merging  criterion  (steepest  descent)  to 
arrive  at  a  minimum  of  E  . 


4.  PROGRAM  STRUCTURE  GENERAL  DESCRIPTION 


The  RG  method  arrives  at  a  segmented  image  by  merging  adjacent  regions  with  a 
positive  merging  criterion  (Equation  6)  until  a  minimum  of  the  simplified  M&Sh 

functional  (Equation  2)  is  reached.  To  evaluate  Equation  6,  we  need  the  coefficient 
vectors  At  and  Aj,  the  two  matrices  M,  and  My,  and  the  length  of  the  common 
boundary  i{d{Ol,Oj)') .  After  merging  two  regions  O,  and  0J}  we  also  need  the  two 
vectors  and  Vj  in  order  to  compute  the  new  coefficient  vector 

A..  =(Mj  +  Mj)~' (Vj  +  Vj)  for  the  approximation  to  g  on  O,  uO,. 

All  the  versions  of  the  RG  algorithm  that  we  have  implemented  keep  track  of  several 
lists  of  objects  that  are  associated  with  each  region.  Each  region  k  is  associated  with  five 
lists  and  three  objects.  The  three  objects  are  the  coefficient  vector  Ak ,  the  vector  Vk ,  and 
the  matrix  Mk  associated  with  region  k. 

The  five  lists  associated  with  region  k  are 

1 .  A  list  of  neighbors  for  region  k 
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2.  A  list  with  the  lengths  of  the  common  boundaries  between  region  k  and  each  of  its 
neighbors 

3.  A  list  with  the  values  of  the  merging  criterion  between  region  k  and  each  of  its 
neighbors 

4.  A  list  with  the  boundary  pixel  labels  between  region  k  and  each  of  its  neighbors 

5.  A  list  of  the  pixels  in  region  k 

These  lists  and  objects  must  be  updated  every  time  two  regions  merge.  Updating  lists 
1 , 4,  and  5  involves  taking  the  union  of  two  sets.  The  pixels  in  the  new  region  formed  are 
simply  the  union  of  the  pixels  in  the  two  merged  regions.  This  union  is  implemented  by 
the  function  merge_pix  and  is  discussed  in  more  detail  in  subsequent  paragraphs.  The  list 
of  boundary  pixels  of  the  new  region  is  the  union  of  the  boundary  pixels  of  the  two 
merged  regions  after  deleting,  from  both  lists,  the  pixels  belonging  to  the  common 
boundary.  This  operation  is  implemented  by  the  function  New_Boundary.  When  two 
regions  merge,  say  region  A  and  region  B,  the  list  of  neighbors  for  the  new  region  is  the 
union  of  the  two  lists  of  neighbors  (list  A  and  list  B)  after  deleting  the  following: 

1 .  B  from  list  A 

2.  A  from  list  B 

3.  All  common  neighbors  from  list  B 

The  operations  1,  2,  and  3  for  updating  the  list  of  neighbors  for  the  new  region  are 
implemented  by  the  function  mrg_N_B  and  the  functions  that  mrgJNMB  calls.  The  new 
region  keeps  the  label  A  and  the  new  list  is  list  A.  Region  B  is  said  to  have  been  absorbed 
by  region  A. 

Every  neighbor  of  a  region  B  lists  B  as  one  of  its  neighbors.  When  a  region  B  gets 
absorbed  by  region  A,  B  no  longer  exists  and  must  be  deleted  from  the  lists  of  every  one 
of  its  neighbors  and  substituted  by  A,  unless  the  neighbor  is  a  neighbor  of  both  A  and  B, 
in  which  case  we  delete  only  B  from  its  list.  This  operation  is  performed  by  the  function 
replace_B_by_A.  This  function  also  updates  the  merging  criterion  because  the  area  of 
the  neighbor  A  now  includes  the  area  of  B  as  well.  As  the  different  lists  of  neighbors  are 
updated,  so  are  the  merging  criteria.  The  common  neighbors  are  dealt  with  by  mrg  N  B. 
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The  length  of  the  boundary  between  two  regions  changes  only  when  the  two  regions 
A  and  B  being  merged  have  a  common  neighbor  C.  In  this  case,  if  A  absorbs  B,  then  the 
length  of  the  boundary  between  the  new  A  and  C  is  simply  the  sum  of  the  lengths  of 
boundary  between  C  and  the  old  A  and  between  C  and  B.  Updating  the  boundary  lengths 
is  also  taken  care  of  by  the  function  mrg_N_B. 

The  coefficient  vectors  A  and  V,  the  matrix  M  and  the  preceding  lists  1  through  4  are 
also  updated  in  mrg_N_B  and  the  functions  that  mrg_N_B  calls. 

Figure  4  shows  the  structure  of  the  PC-RG  program  for  two-dimensional  imagery. 
The  program,  called  Tile_Max_Reg_Grow_B,  is  divided  into  four  parts. 

In  Part  I,  the  five  lists  and  the  three  objects  associated  with  each  initial  region  are 
defined  for  every  region  in  the  initial  segmentation.  In  the  PC  case  the  initial 
segmentation  is  chosen  to  be  the  trivial  segmentation  where  each  region  consists  of  a 
single  pixel.  Thus,  for  an  nxm-image,  there  are  nm  initial  regions  in  the  PC  case.  Each 
region  is  labeled  with  an  integer  k;  1  <  k  <  nm.  The  function  initialize  defines  the  initial 
pixels  in  each  region.  The  function  NGB_PNTR_init_B  defines  the  initial  lists  of 
neighbors,  common  boundary  lengths,  merging  criteria,  boundary  pixels,  and  the  initial 
Ak ,  Vk ,  and  Mk  for  each  region  k  e  {1,2,  •  •  •  ,nm) . 

In  Part  II  of  Tile_Max_Reg_Grow_B,  the  centers  of  the  first  four  levels  of  tiles  are 
defined  by  the  function  tile  and  the  two  functions  that  it  calls.  The  function 
sweep_reg_Tile_B  sweeps  through  the  centers  of  the  tiles  in  order  to  merge  them  with 
their  neighbors  according  to  the  merging  criterion  up  to  the  highest  level  of  tiles  chosen. 

Part  III  finishes  the  merging  of  regions  by  selecting  pairs  of  regions  with  maximal 
merging  criterion  (steepest  descent)  for  merging  until  all  the  remaining  criteria  lie  below 
the  chosen  non-negative  threshold.  The  resulting  sets  of  regions,  boundaries,  and 
coefficients  constitute  the  segmentation  of  the  original  image. 

In  Part  IV,  the  piecewise  approximation  to  the  original  image  is  computed  and  the 
information  is  assembled  properly  for  display.  It  is  possible  to  display,  for  example, 
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individual  regions  or  individual  boundaries  if  desired,  as  well  as  the  complete  set  of 
boundaries  and  the  complete  piecewise  approximation  to  the  image. 

We  note  here  that  Part  II  can  be  skipped  if  one  prefers  a  pure  steepest-descent 
procedure  for  minimizing  Equation  2.  This  may  lead  to  a  different  local  minimum  for 
Equation  2.  Part  II  was  introduced  to  speed  up  the  algorithm. 

This  finishes  the  general  description  of  the  structure  of  our  programs.  Next,  we  go 
into  a  detailed  description  of  each  function  involved. 
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Part  I 


Part  II 


Part  III 


Part  IV 


TiIe_Max_Reg_Grow_B 


initialize 

NGB_  PNTR_  init_  B 


Im_  VecM 
fetch  ADD  NGB 


tile 


sweep_reg_Tile_B 


Sweep_  array  _  New 
Sweep_  array_4 


fetch_ADD_NGB 
mrg_N_B 
merge_  pix 


Elim_  Pnt 
Criterion 

Elim_  and_  Update 
close_  ranks 
New_  Boundary 


Criterion 

^  Update^  criterion 
replace_  B_  by_  A 

{close_  ranks 


f max  criterion 

Elim_  Pnt 
Criterion 

Elim_and_  Update 
close_  ranks 

New_  Boundary  —>  close_  ranks 

merge_  pix 
cls_  rnks 
max  criterion 


I  mrg_  N_B  < 


{Criterion 
Update_  criterion 
replace_B_by_A 


DispIay_Reg_Pix 
Vect_to_Img_l 
Display^  ith_  reg 
Vect_  to_Img_l 
Dispiay_  Boundary 
Vect_to_Img_l 
Display_  ith_  Bndry 
Vect_  to_  Img_J 


FIGURE  4.  The  Piecewise  Constant  Matlab  Program. 
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5.  MATLAB  IMPLEMENTATION  OF  PC-RG  METHOD 


5.1  PART  I.  INITIALIZING 


In  Part  I  we  describe  the  functions  initialize,  NGB_PNTR_init_B,  Im_Vect_l,  and 
fetch  ADD  NGB. 


initialize 


I  NGB_  PNTR_  init_  B 


[  Im_Vect_l 
jfetch_ADD_NGB 


5.1.1  Function  initialize 

The  input  to  this  function  is  nm.  The  outputs  are  the  six  arrays:  Start,  End,  Reg_Pix, 
P_S,  StrtJS,  and  End_S. 


This  simple  function  initializes  six  arrays.  To  understand  the  role  of  these  arrays,  we 
should  explain  here  how  we  have  implemented  the  operation  of  union  of  two  disjoint  sets 
in  Matlab. 

Suppose  we  have  four  disjoint  sets  labeled  A,  B,  C,  D,  and  suppose  the  labels  are 
positive  integers,  say  A  =  l,  B  =  2,  C  =  3,  and  D  -  4.  Suppose  the  elements  of  each  set 
are  numbers  also,  say  A  =  {a,,a2,a3} ,  B  =  {bx,b2,b3} ,  C=  {cl,c2,c3,cA} ,  and 
D  =  {d\,d2  ,d3,dA }  .  To  implement  the  possible  unions  of  these  four  sets  (e.g.  AuC  or 
Byj  D)  we  use  four  arrays.  The  first  array  contains  the  elements  of  all  four  sets;  call  it 
Elements'. 


Elements  =  [a],a2,a3,b],b2,b3,c],c2,c3,cA,dl,d2,d3,dA]. 
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The  second  array  has  four  elements  (the  number  of  sets)  and  indicates  where,  in  the 
array  Elements,  the  members  of  each  set  start;  call  it  Start: 

Start  =  [  1,4,  7,  11  ]. 

So,  the  elements  of  set  A  start  at  Start(  A  )  =  Start(l)  =  1,  the  elements  of  set  B  start 
at  Start(  B  )  =  Start(2 )  =  4,  etc.  Thus,  the  elements  of  Start  are  pointers  to  the  beginning 
of  the  list  of  elements  of  each  set  contained  in  the  array  Elements. 

The  third  array  indicates  where,  in  Elements,  the  members  of  each  set  end;  call  it  End: 

End  =[3,6, 10,  14]. 

So,  the  elements  of  the  set  C  =  3  start  at  Start(C)  =  7  and  end  at  End(C )  =  10. 

The  fourth  array  is  a  pointer  to  the  next  element;  call  it  Pointer.  In  this  example,  we 
have 


Pointer  =  [  2,  3,  0,  5,  6,  0,  8,  9,  10,  0,  12,  13,  14,  0  ]. 

The  elements  of  a  set  can  be  retrieved  from  the  array  Elements  using  the  two  arrays 
Start  and  Pointer.  For  example,  the  elements  of  set  B  =  2  are 

Elements(Start(2 ))  =  Elements^ )  =  bx 

Elements(Pointer(Start(2 )))  =  Elements{Pointer{ 4))  =  Elements (5)  =  b2 
Element  s{Pointer{Pointer{Start(2))))  =  Elements(Pointer(5))  =  b3 

Note  that  Pointer{6)  =  0  and  the  zero  indicates  the  end  of  the  list  for  set  B .  Thus,  the  idea 
is  to  iterate  along  the  array  Pointer,  starting  at  Start(B),  until  the  value  zero  is  reached: 
Start{  B ),  Point er(Start(  B  )),  Pointer(Pointer(Start(  B  ))), ... ,  0.  Also  note  that  for  any  set 
X,  Start(X)  is  not  only  the  index  of  Elements  corresponding  to  the  first  element  of  X,  but 
is  also  the  index  for  starting  the  iteration  along  the  array  Pointer  to  retrieve  the  elements 
of  the  set  X. 
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In  general,  the  arrays  Start,  End,  and  Pointer  are  defined  so  that  Start(X)  is  the  index 
of  the  first  element  of  X,  Point er(Start{Xj)  is  the  index  of  the  second  element  of  X, 
Pointer{Pointer{Sicirt{X)))  is  the  third,  and  so  on, 

Pointer(End(Xj)  -  0  for  all  sets  X  (5.1) 

and  End(X)  is  the  index  of  the  last  element  of  X  When  a  set  X  disappears,  we  set  StartiX) 
=  0  to  indicate  that  it  no  longer  plays  a  role.  One  can  go  through  the  elements  of  a  set  or  a 
list  of  elements  using  the  following  Matlab  while  loop. 


x=Start(C); 
while  x  ~=  0 


Elements(x) 

x=Pointer(x); 


end 


%  The  elements  of  set  C 

%  Pointer  to  the  first  element  in  the  list  C  or  set  C 
%  If  Start  (C)  is  not  zero,  the  set  C  is  still  viable 
%  Ifx  is  not  zero,  there  is  another  element  in  the  set  C 
%  An  element  of  set  C 

%  Advance  the  pointer  to  the  index  of  the  next  element 
%  When  x  -  0  the  list  has  ended 


The  union  of  two  sets,  that  is,  an  assignment  of  the  form  A  =  AuC ,  for  example, 
can  be  implemented  by  concatenating  the  list  for  set  A  with  the  list  for  set  C ,  as  follows. 

1 .  Po  inter (End(  A  ))  =  Start(  C  ) 

2.  End(A)  =  End(C) 

3.  Start(C)  =  0 


Statement  1  concatenates  the  two  lists.  Before  Statement  1,  Pointer{End{A))  =  0 
indicates  the  end  of  list  A .  By  setting  Pointer(End(A ))  =  Start(C),  the  list  continues 
through  the  elements  in  the  list  C .  Statement  2  moves  the  end  of  list  A  to  the  end  of  list 
C.  St  art  { C )  =  0  in  Statement  3,  is  for  bookkeeping  purposes;  it  indicates  that  the  set  C 
has  been  merged  with  another  set  and  no  longer  exists. 
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In  the  two-dimensional  setting,  the  image  is  an  array  of  dimensions  nxm.  The 
integers  n  and  m  denote  the  dimensions  of  the  image  throughout  the  paper  and  the 
programs. 

As  mentioned  in  the  preceding  paragraphs,  the  function  initialize  initializes  six 
arrays.  The  first  three  arrays  ( Reg_Pix ,  Start,  and  End)  are  used  to  keep  track  of  the  pixels 
belonging  to  each  region  and  to  implement  the  unions  of  regions  as  the  merging  of 
regions  proceeds.  Since  the  image  is  nxm,  there  are  nm  pixels.  The  pixels  are  labeled  1 
through  nm  (Figure  5),  starting  with  the  first  row  of  pixels.  Pixel  (i,  j)  receives  the  label 
k  =  (i  - 1  )m  +  j. 


1 

2 

3 

m 

m  + 1 

m  +  2 

m  +  3 

•••  2m 

2m  + 1 

2m  +  2 

2m +  3 

•  ••  3m 

(n  - 1  )m  + 1 

(n  - 1  )m  +  2 

(n  - 1  )m  +  3 

•••  nm 

FIGURE  5.  Pixel  Labels. 


In  the  PC  case,  the  initial  segmentation  is  the  trivial  segmentation.  Each  region 
consists  of  one  pixel;  hence,  there  are  nm  regions.  The  regions  are  labeled  1  through  nm 
as  the  pixels,  so  that  Region(i)  =  {/},/=  1,  2, ... ,  nm. 

The  arrays  Start  and  End  in  initialize,  as  discussed  previously,  indicate  where  the 
pixels  of  each  region  start  and  end.  Since  each  region  has  one  element, 

Start  =  End=  [1, 2,  3, . . . ,  nm  ]. 

The  array  RegPix,  the  pointer  associated  with  the  region  labels  and  pixel  labels,  is  a 
1  x  nm  array  and,  according  to  Equation  5.1  and  the  definition  of  End,  Reg_Pix  =  [0,  0, .  . 
.  ,  0].  Moreover,  the  array  RegPix  also  plays  the  role  of  the  array  Elements.  In  this 
setting  the  array  Elements  is  [1,  2,  3,  ...  ,  nm ],  so  it  is  not  necessary  to  explicitly  define 
it.  To  keep  track  of  the  pixels  in  each  region,  we  need  only  Start,  End,  and  Reg  Pix. 
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The  other  three  arrays  defined  in  the  function  initialize  are  called  P_S,  Strt_S ,  and 
EndJS.  These  are  used  to  keep  track  of  the  non-zero  entries  in  the  array  Start.  A  zero 
entry  in  Start  indicates  that  the  region  corresponding  to  that  entry  has  been  absorbed  by 
(merged  with)  another  region  and  will  no  longer  be  considered.  Thus,  by  keeping  track  of 
the  non-zero  entries  in  Start  one  can  search  efficiently  through  the  remaining  regions.  The 
array  P_S  is  the  pointer,  Strt_S  and  EndJS  are  the  start-  and  end-arrays  associated  with 
the  pointer  PS.  Thus,  initially  StrtJS  =  1,  EndJS  =  nm,  and  P  S  =  [2,  3,  4,  ... ,  nm,  0]. 
The  function  cls_rnks  updates  these  three  arrays  so  that  P  S  will  skip  over  the  zero 
entries  in  Start. 


5.1.2  Function  NGB_PNTR_init_B 

The  function  NGB_PNTR_init_B,  one  of  the  main  functions  in  our  implementation 
of  the  RG  method,  initializes  1 1  arrays: 

1 .  The  array  M,  a  1  x  ( nm )  array 

2.  The  array  V,  a  c  x  (nm)  array 

3.  The  array  C,  a  c  x(nm)  array,  where  c  is  the  number  of  channels 

4.  The  array  Reg_N,  a  3  x  (4nm)  array 

5.  The  array  Start _N,  a  1  x  (nm)  array 

6.  The  array  EndJS,  a  1  x  (nm)  array 

7.  The  array  PlrJJ,  a  1  x  (4 nm)  array 

8.  The  array  Reg_B,  a  1  x  (4 nm)  array 

9.  The  array  Start  JB,  a  1  x  (nm)  array 

10.  The  array  EndJB,  a  1  x  (nm)  array 

1 1 .  The  array  PrtJB,  a  1  x  (4 nm)  array 

c 

With  each  region  k,  a  function  uk  =  ^£jakifi  approximates  the  image  g  on  Ok.  The 

i=i 

coefficient  vector  Ak  =  [ak]  ak2  ■■■ake]r  e  Re  satisfies  MkAk  =  Vk  (see  Proposition  1). 
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In  the  PC  setting,  the  matrix  Mk  -  /u(Ok )  •  / ,  where  ju(Ok)  denotes  the  area  of 
region  Ok .  The  1  x(nm)  array  M  contains  the  areas  of  all  the  initial  regions  which  are  all 
equal  to  one. 

In  the  discrete  case,  when  Ok  is  one  pixel,  say  pixel  k,  the  vector  Vk  =  ^F1  g  d/u 

o„ 

becomes  Vk  =  ^FT(i)g(i)  =  =  g(k)  and  since  Mk=I ,  Ak  =Vk=  g(k),  for  all 

ieOt  i€()k 

k.  Here,  the  cxnm  array  C  contains  all  the  initial  coefficients  Ak  =  g(k )  (1  <  k  <  nm), 
where  c  is  the  number  of  channels.  Thus,  C  -  V  and  V  is  a  c  x  nm  array  with  V(:,k)  = 
g(pixel  k) . 

The  function  Img_Vect_l  transforms  a  1 -channel,  (n  x  m) -image  into  a  (1  x  nm)  array 
V.  For  multichannel  data  we  use  Img_Vect_l  on  each  channel  to  transform  a  ochannel, 
(n  x  m)- image  into  a  (c  x  nm)  array  V as  follows. 

For  multichannel  data  use 
for  i=l:ch 

V(i, :)  =Img_Vect_l  (n,  m,  CHANNEL(i)); 

end 

The  arrays  Start_N,  End_N,  and  Ptr_N  are,  respectively,  the  start,  end,  and  pointer 
arrays  associated  with  the  3  x  (4nm)  array  Reg_N,  which  contains  the  neighbors, 
boundary  lengths,  and  merging  criteria  for  each  region  k.  This  array  plays  the  role  of  the 
array  Elements  discussed  as  an  example  in  the  description  of  the  function  initialize. 

Every  pixel  k  (except  the  pixels  at  the  edges  of  the  image)  has  four  neighboring 
pixels:  an  up,  down,  left,  and  right  neighbor  (see  Section  3).  Each  neighbor  has  a  label  in 
the  set  {1,2,3 ,---,nm) .  The  3  x  4 nm  array  Reg_N  contains,  in  the  first  row,  the  labels  of 
the  neighbors  of  each  pixel.  The  label  0  is  given  to  neighbors  that  do  not  exist.  The 
neighbors  are  labeled  counterclockwise  starting  with  the  upper  neighbor:  up,  left,  down, 
and  right.  So,  for  example  (see  Figure  5),  the  neighbors  of  pixel  1  are  0,  0,  m+ 1,  2.  The 
neighbors  of  pixel  m+2  are  2,  m+ 1,  2m+2,  m+ 3.  The  neighbors  are  listed,  four  per  pixel, 
along  the  first  row  of  Reg_N  in  the  same  order  as  the  pixel  labels. 
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On  the  second  row  of  Reg_N  and  in  the  same  order  are  listed  the  lengths  of  the 
common  boundaries  between  a  region  and  its  four  neighbors  multiplied  by  lambda  (see 
Equation  3).  Since  all  the  regions  consist  of  one  pixel  initially,  all  the  lengths  equal  one. 
Along  the  third  row  of  Reg_N,  and  in  the  same  order  as  the  neighbors,  are  listed  the 
values  of  the  merging  criteria  between  a  region  k  and  its  four  neighbors.  The  initial 
merging  criteria  in  the  PC-setting  are 

lambda  ~^[C(k)  -  C(neighbor)]T[C(k )  -  C(neighbor)] ,  (5.2) 

where  C(k)  are  column  vectors  of  coefficients  (corresponding  to  the  coefficient  vectors 
A i<  of  Equations  3  and  6)  of  dimension  c  (number  of  channels)  associated  with  region  k. 
Since  the  Expression  5.2  involves  an  inner  product  of  c-dimension  vectors,  the  Matlab 
statement  that  implements  it  is  already  multichannel.  However,  one  must  use  C(:,k )  and 
C(:,  neighbor).  Here,  we  call  the  function  fetch_ADD_NGB  to  fetch  the  neighbors  of 
each  region  k.  It  is  interesting  that  the  function  fetch_ADD_NGB  uses  Start_N,  Ptr_N, 
and  Reg_N  as  inputs,  and  these  arrays  are  being  defined  here  in  the  calling  function.  It 
works.  The  parts  that  fetch_ADD_NGB  needs  are  defined  before  it  is  called. 

Since  each  region  has  four  neighbors,  the  list  of  neighbors  for  region  k  ends  at 
Reg_N(l,  4k)  and  starts  at  Reg_N(l,  4k-3)  (1  <  k  <  nm).  Thus,  End_N(k)  =  4k  and 
Start _N(k)  =  4k- 3  (1  <k  <nm).  The  pointer  Ptr_N  is  defined  as 

Ptr_N=  [  2,  3,  4,  0,  6,  7,  8,  0,  10,  1 1,  12,  0, ... ,  4 nm  -2, 4 nm  -1, 4 nm  ,  0]. 

The  arrays  StartJB,  EndJB,  and  Ptr_B  are,  respectively,  the  start,  end,  and  pointer 
arrays  associated  with  the  1  x  (4 nm)  array  Reg_B,  which  contains  the  labels  of  the 
boundaries  associated  with  each  region  k.  This  array  plays  the  role  of  the  array  Elements 
discussed  as  an  example  in  the  description  of  the  function  initialize. 

Between  any  two  adjacent  pixels  we  consider  the  existence  of  a  virtual  boundary 
(virtual  in  the  sense  that  there  is  no  pixel  corresponding  to  the  boundary)  and  divide  the 
boundaries  into  two  classes:  vertical  and  horizontal.  The  vertical  boundaries  are  labeled  1 
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through  nm.  The  horizontal  boundaries  are  labeled  nm  +  1  through  2 nm.  The  number  of 
boundaries  is  fewer  than  the  labels  because  the  boundaries  of  the  image  are  not 
considered  boundaries  of  pixels.  Thus,  not  all  the  labels  are  used. 

Except  for  the  pixels  at  the  edges  of  the  image,  every  pixel  k  has  four  boundaries:  an 
upper  boundary  labeled  (nm+k-m),  a  left  boundary  labeled  (k-1 ),  a  down  boundary 
labeled  (nm+k),  and  a  right  boundary  labeled  k.  The  nonexisting  boundaries  are  labeled  0. 
We  note  here  that  one  can  label  the  horizontal  virtual  boundaries  1  through  nm  as  well. 
Different  labels  are  not  necessary  (see  the  remark  in  Section  5.4.4). 

Reg_B  is  a  lx(4 nm)  array  with  the  labels  of  the  four  boundaries  of  each  pixel  in  the 
same  order  as  the  neighbors.  Thus,  Start_B  =  Start_N,  End_B  =  End_N,  and  Ptr_B  = 
Ptr_N. 


5.1.3  Function  Img  Vect  l 

This  function  transforms  a  single-channel  nxm  image  into  a  (1  x  nm)  vector  V: 

V(l ,(/ - 1  )m  +  j )  =  Image(/,y) ,  1  <i  <n,  1  <  j  <  m. 

In  the  multichannel  setting,  this  function  must  be  modified  to  produce  a  (c  x  nm) 
vector  V,  where  c  is  the  number  of  channels  or  else  the  function  is  used  on  each  channel 
one  at  a  time. 


5.1.4  Function  fetch_ADD_NGB 

This  function  retrieves  the  list  of  addresses  of  the  neighbors  of  a  region  A  and  the  list 
of  neighbors,  boundary  lengths,  and  merging  criteria  using  the  arrays  Start _N,  Ptr_N,  and 
Reg_N.  The  inputs  to  the  function  are  A,  Start _N,  Ptr_N,  and  Reg_N.  The  function  returns 
four  arrays:  an  array  Addresses  with  the  addresses  of  the  neighbors  of  region  A,  an  array 
Neighbs  with  the  labels  of  the  non  zero  neighbors  of  A,  and  the  arrays  Lengths  and  Criter 
with  the  lengths  of  boundaries  and  merging  criteria,  respectively. 
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5.2  PART  II.  TILING 


tile 


sweep_  reg_  Tile_  B 


S\veep_  array  _  New 
Sweep_  array_4 


fetch_ADD_NGB 


mrg_N_  B 


[  merge_  pix 


Elim_Pnt 

Criterion 

Eiim_and__  Update 
close  ranks 


New_  Boundary  ->  {close_  ranks 


Criterion 
Update^  criterion 
replace_B_by_A 


Part  II  of  the  PC-RG  program  in  two  dimensions  implements  the  tiling  described  in 
Section  3,  which  was  introduced  to  speed  up  the  algorithm.  Here  we  describe  the 
functions  tile,  Swcep_array_New,  Sweep_array_4,  and  sweep_reg_Tile_B.  The 
function  fetch_ADD_NGB  was  described  in  Part  I,  and  the  rest  of  the  functions  are 
described  in  Part  III. 


5.2.1  Functions  tile,  Sweep_array_New,  and  Sweep_array_4 

The  function  tile  simply  calls  the  two  functions  Sweep_array_New  and 
Sweep_array_4  to  create  two  arrays,  77  and  T2.  These  arrays  contain  the  centers  of  the 
levels  1,  2,  3,  and  4  tiles  (see  Section  3).  This  function  has  three  inputs:  the  dimensions  of 
the  image  n  and  m  and  level  e  (1,  2,  3,  4},  which  is  used  to  select  the  highest  level  of 
tiles  desired.  If  level  <  2,  then  only  T1  is  created  by  Sweep_array_New.  If  level  >  2,  then 
Sweep_array_4  creates  T2.  The  centers  of  the  first  two-level  tiles  are  defined  in  T1  (the 
first  level  in  row  1  and  the  second  level  in  row  2).  The  centers  of  the  third-  and  fourth- 
level  tiles  are  in  T2  (the  third  level  in  row  1  and  the  forth  level  in  row  2).  The  outputs  of 
this  function  are  T1  and  T2. 
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5.2.2  Function  sweep_reg_Tile_B 


The  inputs  to  this  function  are  n,  m,  M,  V,  C,  T,  Tl,  T2,  level,  S_N,  E_N,  P_N,  R_N,  S, 
E,  RP,  S_B,  E_B,  PJB,  and  R_B. 

This  function,  which  forms  the  tiles,  goes  through  the  centers  of  the  various-level  tiles 
in  Tl  and  T2  according  to  the  level  of  tiling  selected  (input  level)  and  merges  the  center 
of  a  tile  with  its  adjacent  neighbors  (to  form  the  tile),  if  the  merging  criterion  exceeds  the 
chosen  threshold  T.  The  steps  are  as  follows. 

A  center  is  selected  and  labeled  A: 

A  =  T(l,k)  for  level  1,  k  =  1, 2,  3, . . .,  or 

A  =  T(2,k)  for  level  2,  k  =  1, 2, 3, . . .,  or 

A  =  T(3,k)  for  level  3,  k  =  1, 2,  3, . . .,  or 

A  =  T(4,k)  for  level  4,  k=  1,2,  3, . . .. 

If  a  ^  0,  then  sweep_reg_Tile_B  calls  the  function  fetch_ADD_NGB  (described  in 
Section  5.1)  to  get  the  neighbors  of  A.  Then  it  loops  through  the  list  of  neighbors, 
checking  if  the  merging  criterion  exceeds  the  chosen  threshold  T  for  merging,  and  if  it 
does,  it  calls  the  function  mrg_N_B  to  update  the  lists  and  objects  associated  with  A  and 
its  neighbors.  Then  it  calls  the  function  merge_pix  to  implement  the  union  of  the  pixels 
of  A  and  the  neighbor  being  absorbed.  It  returns  the  updated  arrays  M,  V,  C,  S_N,  E_N, 
P_N,  R_N,  S,  E,  R_P,  S_B,  E_B,  PJB,  R_B. 

The  function  mrg_N_B  and  the  functions  it  calls  and  the  function  merge_pix  are 
discussed  in  Part  III. 
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5.3  PART  III.  STEEPEST  DESCENT 


Part  III  implements  the  Steepest  Descent  procedure  for  minimizing  Equation  2. 


max_  criterion 
mrg_N_B 


Elim_  Pnt 
Criterion 


<  merge_pix 


-H 


Elim_  and_  Update 


cls_  mks 
max_  criterion 


close_  ranks 
New_  Boundary 


Criterion 
Update_  criterion 
replace_  B_  by_  A 


—>  {close_  ranks 


5.3.1  Function  maxcriterion 

The  inputs  to  max_criterion  are  Start_S,  Ptr_S,  Start_N ,  Ptr_N,  and  Reg_N.  The 
outputs  are  Region,  Neighbor,  and  MAX_Criter. 

This  function  uses  a  double  while  loop  to  search  through  the  remaining  regions  (using 
Start _S  and  Ptr_S  discussed  in  Section  5.1.1  as  StrtJS  and  P_S)  and  through  each 
remaining  neighbor  of  each  region  (using  Start _N  and  Ptr_N  )  in  order  to  find  the  pair  of 
remaining  neighboring  regions  with  maximal  criterion.  The  function  returns  the  selected 
region-neighbor  pair  and  the  value  of  the  maximal  criterion. 
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5.3.2  Function  mrg_N_B 

This  function  implements  most  of  the  operations  involved  when  two  regions  merge. 
The  inputs  are  A,  B,M,  V,  C,  Start _N,  End_N,  Ptr_N,  Reg_N,  Start _B,  End_B,  Ptr_B ,  and 
Reg_B.  The  two  regions  to  be  merged  are  A  and  B.  The  new  region  formed  keeps  the 
label  A;  that  is,  A  will  absorb  B.  Both  A  and  B  are  integers  between  1  and  nm. 

The  arrays  M,  V,  and  C  contain,  respectively,  the  matrix  Mk ,  the  vector  Vk  and  the 
vector  of  coefficients  Ck  (see  Proposition  1 ;  we  use  Ck  instead  of  Ak )  associated  with 
every  region  k  (1  <  k  <  nm).  In  the  1 -channel  setting,  these  three  objects  are  scalars. 
Thus,  M,  V,  and  C  are  1  x  nm  vectors.  In  Part  1  of  mrg_N_B  the  values  of  these  scalars 
are  computed  for  the  new  region.  The  new  region  keeps  the  label  A.  The  update  of  these 
objects  in  the  single-channel  case  reads  as  follows: 

M(A)  =  M(A)+  M(B),  a  scalar 

V(A)  =  V(A)  +  V(B),  a  scalar 

C(A)  -  V(A)/M(A),  a  scalar 

In  the  multichannel  case  (c  channels),  the  matrices  Mk  are  simply  the  (c  x  c)-identity 
matrix  multiplied  by  the  area  of  the  region.  Therefore,  they  can  be  represented  by  a  scalar 
as  in  the  1  -channel  case.  The  vectors  Vk  and  Ck  have  dimension  c.  Thus,  V  and  C  are  (c 
x  run)  matrices,  where  V( : ,  k)  and  C( : ,  k)  are  column  vectors  corresponding  to  region 
k.  The  update  of  these  objects  in  the  multichannel  case  reads  as  follows: 

M(A )  =  M(A )  +  M(B) ,  a  scalar 

V( : ,  A)  =  V( : ,  A)  +  V( : ,  B),  a  column  vector  of  dimension  c 

C(: ,  A)  =  C(: ,  A  )/M(A),  a  column  vector  of  dimension  c 

Note  that  no  matrix  inversion  is  required.  This  modification  is  the  only  one  needed  in 
mrg  N  B  to  implement  the  multichannel  PC  algorithm. 

In  Parts  2  and  3  of  mrg_N_B,  the  array  RegJJ  is  updated.  As  mentioned  in  Section 
5.1.2,  Reg_N  is  a  3  x  4 nm  array.  The  first  row  of  Reg_N  contains  the  lists  of  neighbors  of 
every  region;  the  second  row  contains  the  lists  of  common  boundary  lengths  between 
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every  region  and  its  neighbors.  The  third  row  contains  the  lists  of  merging  criteria 
between  every  region  and  its  neighbors.  When  talking  about  a  list  of  neighbors,  lengths  or 
criteria  of  a  given  region  A,  we  refer  to  it  as  list(yl)  for  short.  It  will  be  clear  from  the 
context  what  kind  of  list  it  refers  to.  The  length  of  the  common  boundary  between  two 
regions,  A  and  B  say,  is  referred  to  by  length^, B).  Similarly,  the  merging  criterion 
between  two  regions,  A  and  B  say,  is  referred  to  by  criterion^ ,B). 

In  Part  2,  a  double  while  loop  is  used  to  search  through  list(^4)  and  listfB)  for  a 
common  neighbor  CN.  When  a  common  neighbor  CN  is  found,  it  is  deleted  from  list(£) 
to  avoid  duplication  of  CN  in  list(,4  u  B)  when  we  merge  the  two  lists.  Next,  a  call  to  the 
function  Elim_pnt  (described  in  the  following  text)  eliminates  B  from  list(CTV),  for  B  has 
been  absorbed  and  no  longer  exists.  Since  CN  is  a  common  neighbor,  A  does  not  need  to 
be  added  to  list(CTV);  however,  since  CN  is  a  neighbor  common  to  A  and  B,  length^, CN) 
has  increased  and  it  will  have  to  be  updated  in  list(^)  and  list(CvV).  After  the  merger  of  A 
and  B,  we  have 

length^,  CN)  —  length^,  CN)  +  length(£,  CN). 

The  function  Elim_Pnt  also  returns  a  pointer  to  A  in  list  (CN)  in  order  to  update 
length^, CN)  in  list(CTV).  Once  this  pointer  is  available,  length^, CN)  in  list(CA)  is 
updated,  after  updating  length^, CA)  in  list(/f). 

Since  length^,  CN)  has  increased  and  the  area  of  A  has  increased,  the  criterion^,  CN) 
has  changed.  A  call  to  the  function  Criterion  (described  in  a  subsequent  paragraph) 
updates  criterion^,  CN)  in  list(yl).  Then,  using  the  pointer  to  A  in  list  (CN)  previously 
provided  by  Elim_Pnt,  we  update  criterion^, CN)  in  list(CTV).  This  takes  care  of  the 
common  neighbors.  We  suspect  that  it  is  not  necessary  to  update  criterion^, CN)  in 
list(/l)  and  list(CA)  here,  for  it  will  probably  be  taken  care  of  in  Part  3,  where  we  update 
the  criteria  between  A  and  the  rest  of  its  neighbors  that  were  not  common  to  both  A  and  B. 

As  mentioned  in  the  preceding  text,  Part  3  takes  care  of  the  list-updating  for  those 
neighbors  of  A  and  B  that  were  not  common  neighbors.  Part  3  also  updates  the  lists  of  the 
members  of  list (5).  This  is  done  by  a  call  to  the  function  Elim_and_Update  and  is 
explained  subsequently  in  the  description  of  this  function. 


33 


NAWCWD  TP  8525 


Once  all  the  lists  of  have  been  updated  and  the  lists  of  neighbors  of  A  and  B  are  ready 
to  be  merged,  in  Part  4  of  mrg_N_B  the  lists  of  the  neighbors  of  A  and  B  are 
concatenated  as  described  in  Section  5.1.1.  Next,  a  call  to  the  function  close jranks 
adjusts  the  pointer  Ptr_N  to  skip  over  the  addresses  that  correspond  to  regions  in  list(^) 
that  have  been  absorbed  and  adjusts  Start_N  and  End_N.  Finally,  the  boundary  of  the  new 
region  is  adjusted  with  a  call  to  the  function  New_Boundary.  These  last  two  functions 
are  also  subsequently  described. 

5.3.2.1  Function  ElimPnt.  The  inputs  to  the  function  Elim_Pnt  are  C,  B,  A, 
Start _N,  PtrJN,  and  RegJN.  Using  Start_N,  Ptr_N,  RegJN,  and  the  searching  technique 
described  in  Section  5.1.1,  this  function  searches  through  the  list  of  neighbors  of  region  C 
looking  for  regions  B  and  A.  When  it  finds  B,  it  eliminates  B  from  the  list  by  setting  the 
entry  to  zero;  and  when  it  finds  A,  it  saves  the  pointer  to  A  and  returns  it  to  the  calling 
function  (pnt_to_AJnC) . 

53.2.2  Function  Criterion.  The  function  Criterion  has  inputs  A,  C,  pointer,  M, 
Coeff,  and  Reg_N.  The  input  pointer  is  a  pointer  to  region  C  in  the  list  of  neighbors 
list(.4)  that  resides  in  Reg_N.  After  computing  the  new  merging  criterion  between  regions 
A  and  C,  it  places  it  in  the  list  of  criteria  in  list(d)  using  pointer  and  Reg_N(3,:),  namely, 
Reg_N(3, pointer).  Since  Coeff,  the  input  of  coefficients  for  the  piecewise  approximation 
can  be  a  column  vector,  the  function  Criterion  is  already  multichannel.  The  function 
returns  the  updated  array  Reg_N. 

53.23  Function  ElimandUpdate.  Before  merging  regions  A  and  B,  this  function 
updates  list(zf),  the  lists  of  the  neighbors  of  A,  list(B),  and  the  lists  of  the  neighbors  of  B 
as  follows. 

Updating  list(^4)  and  the  lists  of  the  neighbors  of  A, 

1 .  Eliminates  B  as  neighbor  from  the  listed). 

2.  Updates  criterion^, Q  for  every  neighbor  C  in  listfd),  except  neighbor  B,  by 
calling  the  function  Criterion. 

3.  Updates  criterion^, Q  in  list(C)  for  every  neighbor  C  in  list(A),  except  neighbor  B, 
by  calling  the  function  Update_criterion. 
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Updating  list(5)  and  the  lists  of  the  neighbors  of  B, 

1 .  Eliminates  A  as  neighbor  from  the  list(B). 

2.  Updates  criterion^, C)  for  every  neighbor  C  in  list(5),  except  neighbor  A,  by 
calling  the  function  Criterion. 

3.  Replaces  B  by  A  in  list(C)  for  all  the  neighbors  C  of  B  in  list(B)  except  neighbor  A, 
and  updates  criterion^ ,C)  in  list(C)  by  calling  the  function  replace_B_by_A. 

Note  that  if  C  is  a  common  neighbor  of  A  and  B,  then  B  will  have  been  replaced  by  0 
at  the  common  neighbor  stage  (i.e.,  function  Elim_Pnt).  Thus,  A  will  not  be  duplicated  in 
the  list(C). 

5.3.2.3.1  Function  Update  criterion.  The  inputs  to  this  function  are  C,  A,  M,  Coeff, 
Start_N,  Ptr_N,  and  Reg_N.  Using  Start_N,  Ptr_N,  Reg_N,  and  the  searching  technique 
described  in  Section  5.1.1,  the  function  searches  through  the  list  of  neighbors  of  region  C 
looking  for  region  A.  When  it  finds  A  ( Reg_N(l ,  x)  =  A),  it  updates  the  criterion(Cy4)  in 
the  array  Reg_N,  namely,  Reg_N( 3,x)  =  criterion(Cy4),  by  calling  the  function  Criterion. 

5.3.2.3.2  Function  replace_B_by_A.  The  inputs  to  this  function  are  A,  B,  C,  Coeff, 
M,  Start_N,  Ptr_N,  and  Reg_N.  Using  Start_N,  Ptr_N,  Reg_N,  and  the  searching 
technique  described  in  Section  5.1.1,  the  function  searches  through  the  list  of  neighbors 
of  region  C  looking  for  region  B.  When  it  finds  B,  it  replaces  B  by  A  and  updates  the 
criterion(Cy4)  in  list(C)  contained  in  the  array  Reg_N  by  calling  the  function  Criterion. 

5.3.2.4  Function  close  ranks.  The  inputs  to  function  close_ranks  are  A,  Start_N, 
End_N,  Ptr_N,  and  Reg_N.  This  function  modifies  Start_N(A),  End_N(A),  and  the  part  of 
the  array  Ptr_N  that  corresponds  to  the  list  of  neighbors  of  A.  This  is  done  in  order  to  skip 
over  the  addresses  in  the  pointer  Ptr_N  that  correspond  to  the  zeros  in  the  part  of  the 
array  Reg_N(\,  :  )  that  corresponds  to  Iist(/4)  so  that  regions  that  no  longer  exist  (the 
zeros)  are  no  longer  visited  when  searching  the  lists  corresponding  to  region  A.  This 
function  is  called  every  time  a  region  A  merges  with  another  region  B.  When  a  merger 
occurs,  the  list  of  neighbors  of  the  new  region  A  becomes  the  union  of  neighbors  of  A  and 
B  after  deleting  B  from  list(/4),  deleting  A  from  list(J5)  and  any  common  neighbors  from 
list(B)  to  avoid  duplications.  Thus,  at  least  two  new  zero  entries  exist  in  list(^4)  and  must 
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be  skipped  over  by  Ptr_N.  Other  zeros  may  exist  in  list(J)  if  a  neighbor  of  A  has  merged 
with  a  region  that  is  also  a  neighbor  of  A  (see  Elim_Pnt). 

The  first  while  loop  in  close_ranks  modifies  Start Jd(A),  End_N(A),  and  Ptr_N  to 
skip  over  the  zeros  that  might  be  at  the  beginning  of  list(^4)  in  Reg_N(  1,  : ).  The  second 
while  loop  in  close_ranks  modifies  End_N(A),  and  Ptr_N  to  skip  over  the  zeros  that 
might  be  in  list(^4)  after  some  non-zero  entries,  interleaved  between  non-zero  entries 
and/or  at  the  end  of  list(/4). 

S.3.2.5  Function  New_Boundary.  The  inputs  to  New_Boundary  are  A,  B,  Start_B, 
End_B,  Ptr_B,  and  Reg_B.  This  function  forms  the  boundary  of  the  union  of  two  regions 
A  and  B  after  the  merger  of  A  and  B.  There  are  two  steps:  (1)  delete  points  that  are 
common  to  both  boundaries  from  the  two  lists  of  boundary  points  list(^4)  and  listfB),  and 
(2)  concatenate  the  two  lists  and  close  ranks  (skip  over  zero  entries  in  RegJB). 

Step  1 .  Using  a  double  while  loop  the  function  searches  through  list(/l)  and  list(B)  in 
RegJB  for  common  boundary  points,  deleting  them  by  setting  the  corresponding  entries 
in  Reg_B  to  zero. 

Step  2.  The  lists  are  concatenated  using  the  following  three  statements  as  discussed  in 
Section  5.1.1. 

Ptr_B(End_B(A))  =  Start JB(B) 

End_B(A)  =  EndJB(B) 

Start _B(B)  =  0. 

A  call  to  the  function  close_ranks  modifies  Start _B,  End_B  and  the  part  of  the  array 
Plr  B  that  corresponds  to  the  list  of  new  boundary  points  of  A,  to  skip  over  the  addresses 
in  PlrJB  that  correspond  to  zeros  in  the  array  Reg_B  that  resulted  from  the  elimination  of 
common  boundary  points. 
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5.3.3  Function  merge _pix 

The  inputs  to  merge_pix  are  A,  B,  Start,  End,  and  Pointer.  This  function  simply 
concatenates  the  list  of  pixels  of  region  A  with  the  list  of  pixels  of  region  B  as  discussed 
in  Section  5.1.1  and  Step  2  in  the  function  New_Boundary. 


5.3.4  Function  cls  rnks 

The  inputs  to  clsjrnks  are  Strt_S,  End_S,  P_S,  and  Start.  This  function  has  the  same 
structure  and  serves  the  same  purpose  (skip  over  zeros  in  Start)  as  the  function 
close_ranks  previously  described.  It  modifies  StrtJS,  End_S,  and  P_S  in  order  to  skip 
over  the  zeros  in  the  array  Start. 

When  a  merger  occurs,  the  list  of  pixels  of  the  new  region  A  becomes  the  union  of  the 
pixels  of  A  and  the  pixels  of  B  and  the  array  Start  acquires  a  new  zero  at  Start(B).  This 
function  is  called  only  after  every  multiple  of  n  mergers  has  occurred. 
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5.4  PART  IV.  COLLECTING  REGIONS  AND  BOUNDARIES  AND  THE  PC- 
SEGMENTED  IMAGE 


Display_  Reg_  Pix 
Vect_to_Img_l 
Display_ith_reg 
Vect_to_Img_l 
Display_  Boundary 
Vect_to_Img_l 
Display_  ith_  Bndry 
Vect_to_Img_l 


5.4.1  Function  Display_Reg_Pix 

The  inputs  to  the  function  Display_Reg_Pix  are  Start _S,  Ptr_S,  Start,  Reg_Pix,  nm, 
and  C.  The  outputs  are  Regns,  Segmentation,  and  Final  Regions. 

The  arrays  Start_S  and  Ptr_S  are  associated  with  the  labels  of  the  remaining  regions 
(start  and  pointer).  Start _S  is  the  label  of  the  first  remaining  region.  Ptr_S  is  the  pointer  to 
subsequent  regions.  Using  Start _S  and  Ptr_S,  this  function  collects  the  labels  of  all  the 
regions  in  the  final  segmentation  and  puts  them  in  the  array  Final  Regions.  The  size  of 
this  array  is  then  the  number  of  regions  left  in  the  segmentation  obtained. 

The  arrays  Start  and  RegPix  are  associated  with  the  pixels  in  every  remaining  region 
(start  and  pointer).  If  A  =  Start_S,  then  A  is  the  label  of  the  first  region  and  Start(A)  is  the 
first  pixel  in  A.  Reg  Pix  is  the  pointer  to  subsequent  pixels.  Using  Start _S,  PtrjS,  Start 
and  Reg  Pix,  this  function  collects  the  pixels  in  every  remaining  region  and  forms  the 
two  row  vectors  Segmentation  and  Regns  of  dimension  nm.  If  pixel  x  belongs  to  a  region 
A,  then 


Regns(x)  =  A  and  Segmentation (x)  =  C(A)  (single-channel), 
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where  C(A)  is  the  PC  value  of  the  approximation  to  the  original  image  that  is  associated 
with  region  A.  Thus,  after  transforming  the  two  row  vectors  Segmentation  and  Regns  into 
an  nxm-matrix  by  the  function  Vect_to_Img_l  in  Section  5.4.2,  these  two  matrices 
represent  the  remaining  regions  in  the  final  segmentation  (by  label)  and  the  piecewise 
approximation  to  the  original  image,  respectively. 

In  the  multichannel  setting  we  simply  replace  the  last  scalar  assignment  statement 
above  by  a  vector  assignment: 

Segmentation^ ,x)  =  C(:,A)  (multichannel), 

where  C  is  the  cxnm  array  of  piecewise  linear  coefficients.  Then  one  must  transform  the 
cxnm  array  Segmentation  into  c-channel  image  data  by,  for  example,  applying 
Vect_to_Img_l  to  each  of  the  rows  of  Segmentation. 


5.4.2  Function  Vect_to_Img_l 

The  function  Vect_to_Img_l  has  inputs  m,  and  V.  The  dimension  of  the  vector  V  is  a 
multiple  of  m,  say  nm.  The  function  transforms  the  vector  V  into  an  nxm-matrix  called 
Image,  the  output.  The  components  of  the  vector  are  laid  along  the  rows  of  the  matrix. 


5.4.3  Function  Display _ith_r eg 

The  inputs  to  the  function  Display_ith_reg  are  A,  Start,  Reg  Pix,  nm,  and  C.  The 
output  is  the  nm  row  vector  ith_region.  A  is  the  label  of  the  ith-region  obtained  using  the 
array  Final JRegions  provided  by  the  function  Display  Reg  Pix  of  Section  5.4.1.  This 
function  does  for  region  A  what  Display  Reg  Pix  does  for  every  region  in  the  final 
segmentation;  namely,  it  builds  the  nm  row  vector  ithjregion,  which  contains  zeros 
everywhere  except  at  the  locations  corresponding  to  the  pixels  of  region  A.  At  these 
locations  it  will  either  contain  the  label  A  or  the  constant  value  C(A)  of  the  piecewise 
approximation  to  the  original  image  that  is  associated  with  region  A  of  the  segmented 
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image.  The  vector  ithjegion  is  then  transformed  into  an  array  by  the  function 
Vect_to_Img_l  that  can  be  used  to  displays  the  ith-region  by  itself. 

In  the  multichannel  setting,  we  use  the  label  A  or  C(:,A),  the  c- vector  of  coefficients. 
If  the  label  A  is  used,  then  the  output  ithjegion  is  a  row  vector  that  the  function 
Vect  to  Img  1  can  transform  into  an  nxm  image.  If  the  multichannel  coefficient  vector 
C(:,A)  is  used,  then  Vect_to_Img_l  must  be  applied  to  each  row  of  the  cxnm  array 
ithjegion. 


5.4.4  Function  Display  Boundary 

The  inputs  are  Final _Regions,  Start_B,  Ptr_B,  Reg_B,  and  nm.  The  outputs  are  Bndrs 
and  Reg_B.  The  arrays  Start_B  and  Ptr_B  are  the  start  and  pointer  arrays  for  Reg_B, 
which  contains  the  boundary  point  labels.  Reg_B  will  be  modified  by  this  function  as 
described  subsequently. 

There  are  vertical  virtual  boundaries  (see  Section  5.1.2  about  virtual  boundary  points) 
between  consecutive  pixels  on  the  same  row  with  labels  1  through  nm.  To  display  a 
vertical  virtual  boundary  point,  we  use  the  pixel  with  the  same  label.  Thus,  this  is  a  left 
representation  in  the  sense  that  the  boundary  between  two  consecutive  horizontal  pixels  is 
represented  by  the  pixel  on  the  left. 

The  horizontal  virtual  boundaries,  labeled  {nm  +  1)  through  2 nm,  are  the  boundaries 
between  two  consecutive  pixels  on  the  same  column.  To  display  a  horizontal  virtual 
boundary  point,  we  use  the  pixel  with  label  equal  to  the  label  of  the  boundary  point  minus 
nm.  Thus,  this  is  an  up  representation  in  the  sense  that  the  boundary  between  two 
consecutive  vertical  pixels  is  represented  by  the  pixel  on  the  top. 

Remark.  Giving  the  horizontal  virtual  boundaries  a  label  different  from  the  labels  of 
the  vertical  ones  is  not  necessary.  At  the  end  (at  this  point)  we  end  up  using  the  same 
pixel  in  the  representation  of  the  boundary  points  whose  labels  differ  by  nm.  Thus,  a 
different  label  is  not  necessary  to  begin  with  (see  Section  5.1.2.). 
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This  function  builds  the  binary  nm- row  vector  Bndrs  with  a  one  at  every  entry 
corresponding  to  a  pixel  representing  a  boundary  point  according  to  the  preceding 
convention. 

The  array  Final  Regions  provides  the  labels  of  the  regions  in  the  final  segmentation. 
If  A  =  Final  Regions  (i),  then  A  is  the  label  of  the  ith-region.  Start _B (A)  is  the  address  in 
RegJB  of  the  first  boundary  point  of  A;  that  is,  Reg_B(Start_B(A))  is  the  first  boundary 
point  of  A.  Using  a  while  loop  inside  a  for  loop  that  sweeps  through  all  the  final  regions, 
the  function  searches  through  the  array  of  boundary  labels  in  Reg_B.  Only  if  Reg_B(x)  > 
nm,  is  Reg_B  modified  by  setting  Reg_B(x)  =  Reg_B(x)  -  nm.  Then,  Bndrs (Reg_B(x))  is 
set  equal  to  1 . 


5.4.5  Function  Display_ith_Bndry 

The  function  Display_ith_Bndry  builds  the  binary  nm- row  vector  ith_Boundary,  the 
output  of  the  function,  which  contains  a  one  at  every  entry  corresponding  to  a  pixel 
representing  a  boundary  point  of  the  ith-region  in  the  final  segmentation  according  to  the 
convention  in  the  description  of  the  previous  function.  It  can  be  used  to  display  the 
boundary  of  that  region  alone.  The  inputs  to  the  function  are  A,  Start _B,  PtrJB,  Reg_B, 
and  nm. 
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Appendix 

MATLAB  PROGRAM 


function  [Regs,Segmt,Boundaries]  =  Tile_Max_Reg_Grow_B(n,m,level,Image, Threshold, lam) 
%  PIECEWISE  CONSTANT,  ONE  CHANNEL. 

T=Threshold; 

nm=n*m; 

lambda=larn*lam*sqrt(nrn); 

'INITIALIZING'  %  PART  1 . 

[S,E,RP,P_S,Strt_S,End_S]=initialize(nm); 

[M,V,C,S_N,E_N,P_N,S_B,E_B,P_B,R_N,R_B]=NGBJPNTR_init_B(Image,lambda,n,m); 

'TILING'  %  PART  2. 

[Tl,T2]=tiIe(n,m, level); 

[M,V,C,S_N,E_N,P_N,R_N,S,E,RP,S_B,E_B,P_B,R_B]=swp_rg_Tile_B(n,m,M,.... 
V,C,T,T1,T2, level, S_N,E_N,P_N,R_N,S,E,RP,S_B,E_B,P_B,R_B); 

'STEEPEST  DESCENT  %  PART  3. 
[R,N,MAX_Criter]=max_criterion(Strt_S,P_S,S_N,P_N,R_N); 
iteration  =  -1; 

while  MAX_Criter  >  Threshold 
iteration=iteration+ 1 ; 

[M,V,C,S_N,E_N,P_N,R_N,S_B,E_B,P_B,R_B]=mrg_N_B(R,N,M,V,C,... 

S_N,E_N,P_N,R_N,S_B,E_B,P_B,R_B); 

[S,E,RP]=merge_pix(R,N,S,E,RP); 
if  rem(iteration,fix(n))  —  0 
[Strt_S,End_S,P_S]=cIs_rnks(Strt_S,End_S,P_S,S); 

ITERATION=iteration 

end 

[R,N,MAX_Criter]=max_criterion(Strt_S,P_S,S_N,P_N,R_N); 

end 

%  PART  4. 
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'COLLECTING  REGIONS,  BOUNDARIES,  &  THE  PC-SEGMENTED  IMAGE' 

[Regs,Segmt,Final_Regions]=DispIay_Reg_Pix(Strt_S,P_S,S,RP,C); 

Number_of_Regions=length(Final_Regions) 

Regs=Vect_to_Img_l(m,Regs); 

Segmt=Vect_to_Img_l(m,Segmt);  %  SINGLE-CHANNEL 
%for  i=l:length(Final_Regions) 

%  ith_region=DispIay_ith_reg(Final_Regions(i),S,RP,nm,C); 

%  ith_region=Vect_to_Img_l(m,ith_region); 

%  pause 
%end 

[Boundaries,R_B]=Display_Boundaries(Final_Regions,S_B,P_B,R_B,nm); 
Boundaries=Vect_to_Img_l(m, Boundaries); 

%for  i=l:length(Final_Regions) 

%  ith_Boundary=Displ_ith_Bndry(FinaI_Regions(i),S_B,PJB,R_B,nm); 

%  ith_Boundary=Vect_to_Img_l(m,ith_Boundary); 

%  pause 
%end 

subplot(3,l,2),imagesc(Segmt);  axis  image;  %mesh(Segmt) 
subpIot(3,l,3),imagesc(Boundaries);  axis  image;  %mesh(Boundaries) 


PARTI 

INITIALIZING 


function  [Start, End,Reg_Pix, PS, Strt_S,End_S]  =  initialize(nm) 

%  To  initialize  Start,  End,  Reg  Pix,  P_S,  Strt_S,  and  End_S 

Reg_Pix=zeros(l,nm);  %  Pointer  to  pixels  in  region  k.  Pixels  start  at 
Start=[  1  :nm];  %  Start(k)  and  continue  along  Reg_Pix(*)  by 

End=Start;  %  iterating:  Reg_Pix(Start(k)), 

%  Reg_Pix(Reg_Pix(Start(k))),.. .untill  we  reach 
%  a  zero;  i.e.  Reg_Pix(End(k))=0. 

P_S=[[2:nm]  0];  %  Pointer  to  the  non-zero  entries  in  Start. 

Strt_S=l;  %  Beginning  of  non-zero  entries  in  Start. 

End_S=nm;  %  End  of  non-zero  entries  in  Start. 


function  [M, V,C, Start_N,End_N,Ptr_N, Start_B,End_B , .. . 
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Ptr_B,Reg_N,Reg_B]  =  NGB_PNTR_init_B(Image, lambda, n,m) 

%  To  initialize  M,  V,  C,  Start_N,  End_N,  Ptr_N,  Start_B,  End_B,  Ptr_B, 

%  RegN,  and  Reg_B. 

%  Reg_N(l,:)=List  of  neighbors. 

%  Reg_N(2,:)=List  of  common  boundary  lengths. 

%  Reg_N(3,:)=List  of  merging  criteria. 

%  Reg_B=List  of  common  boundary  points. 

%  Start_N=beginnings  of  neighbor  lists. 

%  Start_B=beginnings  of  boundary  points  lists. 

%  End_N=ends  of  neighbor  lists. 

%  End_B=ends  of  boundary  points  lists. 

%  Ptr_N=pointer  to  neighbors  in  the  lists  in  Reg_N(l,:),  a  lx4nm-array. 

%  Ptr_B=pointer  to  boundary  points  in  the  lists  in  Reg_N(4,:). 

%  Initially  Start_N=Start_B,  End_N=End_B,  and  Ptr_N=Ptr_B. 

nm=n*m; 

Q=ones(l,nm); 

M=Q; 

V=Img_Vect_l  (n,m, Image);  %  For  multichannel  data  use: 

%  for  i=l:ch 

%  V(i,:)=Img_Vect_l(n,m,CHANNEL(i)); 

%  end 

C=V; 

End_N=4*[l:nm]; 

Start_N=End_N-3  *Q; 

Ptr_N(End_N)=zeros(  1  ,nm); 

Ptr_N(Start_N)=Start_N+Q; 

Ptr_N(Start_N+Q)=Start_N+2*Q; 

Ptr_N(Start_N+2  *  Q)=Start_N+3  *Q; 

Start_B=Start_N; 
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End_B=End_N; 

Ptr_B=Ptr_N; 

Reg_N=zeros(3,4*nm); 

Reg_B=zeros(l,4*nm); 


for  i=2:n-l 

forj=2:m-l  %  u,  1,  d,  r. 
k=(i-l)*m+j; 
k4=4*k; 

Reg_N(l,k4-3)=k-m;  %  u(k) 

Reg_N(2,k4-3)=l; 

Reg_N(  1  ,k4-2)=k- 1 ;  %  l(k) 

Reg_N(2,k4-2)=l; 

Reg_N(  1  ,k4- 1  )=k+m;  %  d(k) 

Reg_N(2,k4- 1  )=1 ; 
Reg_N(l,k4)=k+l;  %r(k) 
Reg_N(2,k4)=l; 

Reg_B(k4-3)=nm+k-m;  %  u(k) 

Reg_B(k4-2)=k- 1 ;  %  l(k) 

Reg_B(k4-l)=nm+k;  %  d(k) 
Reg_B(k4)=k;  %  r(k) 

end 

%j=m  u,  1,  d. 
k=i*m; 
k4=4*k; 

Reg_N(  1  ,k4-3)=k-m;  %  u(k) 

Reg_N(2,k4-3)=l; 
Reg_N(l,k4-2)=k-l;  %  l(k) 

Reg_N(2,k4-2)=l; 

Reg_N(  1  ,k4- 1  )=k+m;  %  d(k) 

Reg_N(2,k4-l)=l; 

Reg_B(k4-3)=nm+k-m;  %  u(k) 

Reg_B(k4-2)=k- 1 ;  %  l(k) 

Reg_B(k4- 1  )=nm+k;  %  d(k) 

%j=l  u,  d,  r. 
k=(i-l)*m+l; 
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k4=4*k; 

Reg_N(  1  ,k4-3)=k-m;  %  u(k) 

Reg_N(2,k4-3)=l; 

Reg_N(  1  ,k4- 1  )=k+m;  %  d(k) 

Reg_N(2,k4*l)=l; 

Reg_N(  1  ,k4)=k+ 1 ;  %  r(k) 

Reg_N(2,k4)=U 

Reg_B(k4-3)=mn+k-m;  %  u(k) 

Reg_B(k4-l)=nm+k;  %  d(k) 
Reg_B(k4)=k;  %  r(k) 


end 

%  i=l  1,  d,  r. 
for  j=2:m-l 

k=j; 

k4=4*k; 

Reg_N(l,k4-2)=k-l;  %  l(k) 

Reg_N(2,k4-2)=  1 ; 

Reg_N(  1  ,k4- 1  )=k+m;  %  d(k) 

Reg_N(2,k4-l)=l; 
Reg_N(l,k4)=k+l;  %  r(k) 
Reg_N(2,k4)=l; 

Reg_B(k4-2)=k- 1 ;  %  l(k) 

Reg_B(k4-l)=nm+k;  %  d(k) 
Reg_B(k4)=k;  %  r(k) 

%  i=n  u,  1,  r. 
k=(n-l)*m+j; 
k4=4*k; 

Reg_N(  1  ,k4-3)=k-m;  %  u(k) 

Reg_N(2,k4*3)=l; 
Reg_N(l,k4-2)=k-l;  %  l(k) 

Reg_N(2,k4-2)=l ; 

Reg_N(  1  ,k4)=k+ 1 ;  %  r(k) 

Reg_N(2,k4)=l; 

Reg_B(k4-3)=nm+k-m;  %  u(k) 

Reg_B(k4-2)=k-l;  %  l(k) 
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Reg_B(k4)=k;  %  r(k) 

end 

%(ij)=(l,l),k=l,  d,  r. 

Reg_N(l,3)=m+l;  %  d(k) 

Reg_N(2,3)=l; 

Reg_N(l,4)=2;  %r(k) 

Reg_N(2,4)=l; 

Reg_B(3)=nm+ 1 ;  %  d(k) 

Reg_B(4)=l;  %  r(k) 

%  (ij)=(l,m),  k=m,  1,  d. 

k4=4*m; 

Reg_N(  1  ,k4-2)=m- 1 ;  %  l(k) 

Reg_N(2,k4-2)=l; 
Reg_N(l,k4-l)=m+m;  %  d(k) 
Reg_N(2,k4-l)=l; 

Reg_B(k4-2)=m- 1 ;  %  l(k) 

Reg_B(k4-l)=nm+m;  %  d(k) 

%(ij)=(n,l)  u,  r. 

k=(n-l)*m+l; 

k4=4*k; 

Reg_N(l,k4-3)=k-m;  %  u(k) 

Reg_N(2,k4-3)=l; 

Reg_N(  1  ,k4)=k+ 1 ;  %  r(k) 

Reg_N(2,k4)=l; 

Reg_B(k4-3)=nm+k-m;  %  u(k) 

Reg_B(k4)=k;  %  r(k) 

%  (iJ)=Cn,m)  u,  1. 

k=nin; 

k4=4*k; 
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Reg_N(l,k4-3)=k-m;  %  u(k) 

Reg_N(2,k4-3)=l ; 

Reg_N(l,k4-2)=k-l;  %  l(k) 

Reg_N(2,k4-2)=l; 

Reg_B(k4-3)=nm+k-m;  %  u(k) 

Reg_B(k4-2)=k-l;  %  l(k) 

Reg_N(2,:)=lambda*Reg_N(2,:); 

for  k=l  :nm  %  To  compute  the  Merging  Criterion(A,B)  for  all  A=1 ,2,...,nm 
%  and  all  neighbors  B  of  A.  A=k  and  B=Reg_N(l,Ad(i)). 
[Ad,N,L,Crit]=fetch_ADD_NGB(k,Start_N,Ptr_N,Reg_N); 
if  length(Ad)  >=  1 
len=length(Ad); 
for  i=l:len 

M_AB=M(k)+M(Reg_N(l  ,Ad(i)));  %  B=Reg_N(l,Ad(i))  is  a  nghbr  of  A. 

M_i  n  v_ A  B = 1  /M_ A  B ; 

H=M(k)*M_inv_AB*M(Reg_N(  1  ,Ad(i))); 

Reg_N(3,Ad(i))=lambda-(C(k)-C(Reg_N(l)Ad(i)))rH*(C(k)-C(Reg_N(l,Ad(i)))); 

%  FOR  MULTICHANNEL  DATA  USE 

%Reg_N(3,Ad(i))-lambda-(C(:,k)-C(:,Reg_N(l,Ad(i))))'*H*(C(:,k)-C(:,Reg_N(l,Ad(i)))); 

end 

end 

end 


function  V  =  Img_Vect_l(n,m,Image) 

%  This  function  transforms  a  matrix  representation  of  an  image  to  a  vector 
%  representation  of  the  image  by  concatenating  the  rows  of  the  nxm-matrix 
%  to  form  an  nm-row-vector  V. 

%  Jorge  M.  Martin,  NAWC-WPNS  Code  474400D,  China  Lake,  CA.  April  1996. 


for  i=l:n 
i_l=i-l; 
forj=l:m 

V(  1  ,i_l  *m+j)=Image(i,j); 
end 
end 
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function  [Addresses,Neighbs, Lengths, Criter]  =  fetch_ADD_NGB(A,Start_N,Ptr_N,Reg_N) 

%  To  get  the  addresses,  labels,  neighbors,  boundary  lengths, 

%  and  criteria  of  region  A. 

%  Addresses  are  the  addresses  of  the  labels  of  the  neighbors  of  A. 

%  Neighbs  are  the  labels  of  the  neighbors  of  A. 


if  Start_N(A)  —  0 

'A  has  no  neihgbors,  A  has  been  absorbed.' 
return 
end 


x=Start_N(A); 
while  x  ~=  0 

if  Reg_N(l,x)  ~=  0 
Addresses=[x]; 
Neighbs=[Reg_N(  1  ,x)]; 
Lengths=[Reg_N(2,x)]; 
Criter=[Reg_N(3,x)]; 
end 

x=Ptr_N(x); 


%  First  remaining  neighbor  of  A. 

%  The  address  of  the  first  neighbor. 
%  The  first  neighbor. 

%  The  first  bouncary  length. 

%  The  first  merging  criterion. 

%  Next  neighbor. 


PART  II 
TILING 


function  [T1,T2]  =  tile(n,m, level) 

%  To  define  the  centers  of  the  first  four  level  tiles. 
%  level=l,2,3  or  4. 

T 1  =sweep_array_New(n,m); 
if  level  <=  2 
T2=[]; 
return 
end 


A- 8 


NAWCWD  TP  8525 


T2=sweep_array_4(n,m,T  1  (2, :)); 


function  Swp_Arr  =  sweep_array_New(n,m) 

%  To  define  the  centers  of  the  first  two  level  tiles. 
%  Swp_Arr=sweep_array_New(  18,12) 

a=zeros(l,5); 

a(l)=fix(m/5); 

a(2)=fix((m+3)/5); 

a(3)=fix((m+l)/5); 

a(4)=fix((m+4)/5); 

a(5)=fix((m+2)/5); 

b=fix(n/5); 

c=rein(n,5); 

dim=b*sum(a); 
if  c  >=  1 

dim=dim+sum(a(  1  :c)); 
end 

Swp_Arr=zeros(2,dim); 

S  wp_Arr(l ,  1  :a(  1  ))=5  *  [  1  :a(  1 )]; 
aa=a(  1  )+a(2); 

Swp_Arr(l  ,a(l  )+l  :aa)=(m-3)+5*[  1  :a(2)]; 
aaa=aa+a(3); 

Swp_Arr(l  ,aa+l  :aaa)=(2*m-l)+5*[l  :a(3)]; 


if  n  >=  3 

Swp_Arr(2,l:a(3))=Swp_Arr(l,aa+l:aaa);  %  First  block  of  Swp_Arr(2,:). 
end 

aa=aaa; 

aaa=aa+a(4); 

Swp_Arr(l  ,aa+l  :aaa)=(3*m-4)+5*[l  :a(4)]; 
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aa=aaa; 

aaa=aa+a(5); 

Swp_Arr(l  ,aa+l  :aaa)=(4*m-2)+5*[l  :a(5)]; 

if  b  >  1 
for  i=l:b-l 

Swp_Arr(l  ,i*aaa+l  :(i+l)*aaa)=5*m+Swp_Arr(l,(i-l)*aaa+l  :i*aaa); 
end 
end 


if  c  >=  1 
aa=b*sum(a); 
aaa=aa+sum(a(l  :c)); 

Swp_Arr(  1  ,aa+ 1  :aaa)=b*5*m+Swp_Arr(  1 , 1  :aaa-aa); 
end 


%  The  rest  of  Swp_Arr(2,:). 
if  n  >=  8 

d=fix((n-3)/5);  %d>=l. 

for  i=l:d 

Swp_Arr(2,i*a(3)+l:(i+l)*a(3))=(2+5*i)*m-l+5*[l:a(3)];  %  a(3)-Blocks. 
end 
end 


function  Swp  Arr  =  sweep_array_4(N,M,Arr) 

%  To  define  the  centers  of  the  level  3  and  4  tiles. 

m=fix((M+l)/5); 

n=l+fix((N-3)/5); 

Array=Arr(l  :n*m);  %  <====  Not  necessary  ??. 
a=zeros(l,5); 

a(l)=fix(m/5); 

a(2)=fix((m+3)/5); 

a(3)=fix((m+l)/5); 

a(4)=fix((m+4)/5); 
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a(5)=fix((m+2)/5); 

b=fix(n/5); 

c=rem(n,5); 

if  a(l)  =  0 
return 
end 

dim=b*sum(a); 
if  c  >=  1 

dim=dim+sum(a(l:c)); 

end 

Swp_Arr=zeros(2,dim); 

Swp_Arr(  1 , 1  :a(  1  ))=Array(5  *[  1  :a(  1 )]); 
aa=a(  1  )+a(2); 

S  wp_Arr(  1  ,a(  1 )+ 1  :aa)=Array((m-3)+5  *  [  1  :a(2)]); 
aaa=aa+a(3); 

S  wp_Arr(  1  ,aa+ 1 :  aaa)=Array((2  *m- 1  )+5  *  [  1 :  a(3 )] ); 


if  n  >=  3 

Swp_Arr(2,l  :a(3))=Swp_Arr(l,aa+l :aaa);  %  First  block  of  Swp_Arr(2,:). 
end 


aa=aaa; 

aaa=aa+a(4); 

S  wp_Arr(  1  ,aa+ 1  :aaa)=Array((3  *m-4)+5  *  [  1  :a(4)]); 

aa=aaa; 

aaa=aa+a(5); 

Swp_Arr(l,aa+l  :aaa)=Array((4*m-2)+5*[l  :a(5)]); 

if  b  >  1 
for  i=l:b-l 

Swp_Arr(  1  ,i*aaa+l  :(i+ 1  )*aaa)=25  *M+Swp_Arr(  1  ,(i- 1  )*aaa+ 1  :i*aaa); 
end 
end 
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if  c  >=  1 
aa=b*sum(a); 
aaa=aa+sum(a(  1  :c)); 

Swp_Arr(  1  ,aa+ 1  :aaa)=b*25  *M+Swp_Arr(  1 , 1  :aaa-aa); 
end 


%  The  rest  of  Swp_Arr(2,:). 
if  n  >=  8 

d=fix((n-3)/5);  %d>=l. 

for  i=l:d 

S  wp_Arr(2,  i  *  a(3 )+ 1  :(i+l)*a(3))=25*i*M+Swp_Arr(2,l  :a(3));  %  a(3)-Blocks. 

end 
end 


function  [M,V,C,SJ4,EJ4,P_N,R_N,S,E,RP,S_B,EJB,P_B,R_B]  =  swp_rg_TUe_B(n,m,M,.... 
V,C,T,T  1  ,T2, level, S_N,E_N,P_N,R_N,S,E,RP,SJB,E_B,P_B,RJB) 

%  level=l,2  3,  or  4. 

nm=n*m; 

J=1 

if  level  >=  2 
J=2; 
end 

for  j=l  :J 

for  k=  1  :length(T  1(1,:)) 

A=Tl(j,k);  %  Region  A. 

if  A~=0 

if  S_N(A)  ~=  0  %  Region  A  has  not  been  absorbed  yet. 

[Addresses, Neighbs, Lengths, Criter]=fetch_ADD_NGB(A,S_N,P_N,R_N); 
L=length(Neighbs); 
if  L  >=  1 
fori=l:L 
B=Neighbs(i); 

if  Criter(i)  >  T  %  Merge  A  and  B=Neighbs(i). 

[M,V,C,S_N,E_N,P_N,R_N,S_B,E_B,P_B,R_B]=mrg_N_B(A,B,M,V,C,... 
S_N,E_N,P_N,R_N,S_B,E_B,P_B,R_B); 


A-12 


NAWCWD  TP  8525 


[S,E,RP]=merge_pix(A,B,S,E,RP);  %  RP  is  a  pointer, 
end 
end 
end 
end 
end 
end 
end 

if  level  >=  3 
J=l; 

if  level  >=  4 
J=2 
end 

for  j=l  :J 

for  k=l  :length(T2(l,:)) 

A=T2(j,k);  %  Region  A. 

if  A  ~=  0 

if  S_N(A)  ~=  0  %  Region  A  has  not  been  absorbed  yet. 

[A  ddresses,Neighbs, Lengths, Criter]=fetch_ADD_NGB(A,S_N,P_N,R_N); 
L=length(Neighbs); 
if  L  >=  1 
for  i=l:L 
B=Neighbs(i); 

if  Criter(i)  >  T  %  Merge  A  and  B=Neighbs(i). 

[M,V,C,S_N,E_N,P_N,R_N,S_B,E_B,P_B,R_B]=mrg_N_B(A,B,M,V,C,... 

S_N,E_N,P_N,R_N,S_B,E_B,P_B,R_B); 
[S,E,RP]=merge_pix(A,B,S,E,RP);  %  RP  is  a  pointer, 
end 
end 
end 
end 
end 
end 
end 
end 
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PART  III 

STEEPEST  DESCENT 

function  [Region, Neighbor, MAX_Criter]  =  max_criterion(Start_S,Ptr_S,Start_N,... 

Ptr_N,Reg_N) 

%  To  find  the  max  criterion  and  the  region  and  neighbor  corresponding  to  it 


Region=0; 

Neighbor=0; 

MAX_Criter=-l; 

k=Start_S; 
while  k  ~=  0 
x=Start_N(k); 
while  x  ~=  0 

if  Reg_N(l,x)~=  0 
if  MAX_Criter  <  Reg_N(3,x) 
MAX_Criter=Reg_N(3  ,x); 
Neighbor=Reg_N(l  ,x); 
Region=k; 
end 
end 

x=Ptr_N(x); 

end 

k=Ptr_S(k); 

end 


%  First  non-0  region.  Replaces  "for  k=l:nm 

%  First  neignbor  of  region  k. 

%  A  non  zero  neighbor. 

%  Merging  criterion  between  k  &  x. 

%  Merging  criterion  is  larger. 

%  Neigbor  with  largest  citerion  yet. 

%  Region  with  largest  citerion  yet. 

%  Next  neighbor. 

%  Next  region  k  until  k=0. 


function  [M,V,C,Start_N,End_N,Ptr_N,Reg_N,Start_B,End_B,Ptr_B,Reg_B]  =  mrg_N_B(A,... 
B,M,V,C,Start_N,End_N,Ptr_N,Reg_N,Start_B,End_B,PtrJB,Reg_B) 

%  To  update  the  pointer  "Ptr_N"  and  the  arrays  Start_N,  End_N,  and  "Reg_N" 

%  when  the  merger  of  two  regions  A  and  B  occurs. 

%  The  pointer  points  to  "Reg_N"  which  contains  the  labels  of  the  neighbors 
%  of  each  region.  When  regions  A  and  B  are  merged,  the  new  region  formed  is 
%  labeled  region  A.  The  addresses  of  the  neighbors  of  A  start  at 
%  "Start_N(A)"and  end  at  "End_N(A)".  They  are:  Start_N(A), 

%  Ptr_N(Start_N(A)),  Ptr_N(Ptr_N(Start_N(A))), Ptr_N(End_N(A))=0  marks 
%  the  end  of  the  list  for  A. 
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% 

%  C  is  a  common  neighbor  of  A  and  B  which  are  to  be  merged;  A  will  absorb  B 
%  To  eliminate  B  as  a  neighbor  from  the  list(C),  eliminate  C  from  the 
%  list(B)  to  avoid  duplication  of  C  in  the  merged  list(AUB),  and  provide  a 
%  pointer  to  A  in  the  list(C)  to  increment  length(A,C)  in  list(C)  after  the 
%  merging  of  A  and  B. 

% . 

%  PART  1.  PC-2D-lCh.  and  %  Multichannel  Changes. 

%  Update  M(A),  V(A),  and  C(A). 

M(A)=M(A)+M(B);  %  a  scalar.  Single  &  Multichannel. 

M_inv_A=l/M(A);  %  a  scalar.  Single  &  Multichannel. 

V(A)=V(A)+V(B);  %  V(:,A)=V(:,A)+V(:,B);  Multichannel. 

C(A)=M_inv_A*V(A);  %  C(:,A)=M_inv_A*V(:,A);  Multichannel. 

% . 

%  PART  2. 

x=Start_N(A); 
y=Start_N(B); 
while  x  ~=  0 
while  y  ~=  0 

if  Reg_N(l,x)  ~=0&  Reg_N(l,x)  =  Reg_N(l,y) 

%  Found  a  common  neighbor  CN  =  Reg_N(l,x). 
CN=Reg_N(  1  ,x);  %  x  =  pointer  to  CN  in  list(A) 

%  y  =  pointer  to  CN  in  list(B) 

Reg_N(l,y)=0;  %  Delete  CN  from  list(B)  to  avoid 

%  duplication  in  the  merged  list(AUB). 

[Pnt_to_A_inC,Reg_N]=Elim_Pnt(Reg_N(l,x),B)A,Start_N,Ptr_N,Reg_N); 

%  Eliminate  B  from  list(C);  (B  will  be  absorbed 
%  by  A),  and  provide  pointer_to_A_in_C  to  update 
%  length(A,C)  and  Criterion(A,C)  in  list(C). 

RegJN (2 ,x)=Reg_N (2 ,x)+Reg_N (2,y ) ;  %  Update  length(A,C)  in  list(A): 

%  length(A,C)=length(A5C)+length(B,C). 
Reg_N(2,Pnt_to_A_inC)=Reg_N(2,x);  %  Update  length(A,C)  in  list(C): 
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%  length(A,C)=length(A,C)+length(B,C). 
Reg_N=Criterion(A,CN,x,M,C,Reg_N);  %  New  merging  criterion 

%  between  A  and  CN  put  in  list(A). 
Reg_N(3,Pnt_to_A_inC)=Reg_N(3,x);  %  New  merging  criterion 

%  between  A  and  C  put  in  list(C). 

end  %  if 
y=Ptr_N(y); 
end  %  while  y  ~=  0 
x=Ptr_N(x); 
y=Start_N(B); 
end  %  while  x~=  0 

% . 


%  PART  3. 

Reg_N=Elim_and_Update(A,B,C,M,Start_N,Ptr_N,Reg_N); 

%  Update  list(A),  list(B),  and 
%  the  lists  of  the  members  of  list(B) 
%  before  concatenating  the  two  lists. 

% . 


%  PART  4. 

Pn_N(End_N(A))=Start_N(B);  %  To  concatenate  the  lists  for  A  and  B. 

End_N(A)=End_N(B);  %  The  end  of  the  list  for  A  is  now  at  the  end  of  B. 

Start_N(B)=0;  %  Region  B  has  been  absorbed  by  A. 

[Start_N,End_N,Ptr_N]=close_ranks(A,Start_N,End_N,Ptr_N,Reg_N); 

[Start_B,End_B,Ptr_B,Reg_B]=New_Boundary(A,B.Start_B,End_B,Ptr_B,Reg_B); 


function  [Pnt_to_A_inC,Reg_N]  =  Elim_Pnt(C,B,A,Start_N,Ptr_N,Reg_N) 


%  Partially  updating  the  list(C). 

%  C  is  a  common  neighbor  of  A  and  B  which  are  to  be  merged;  A  will  absorb  B. 

%  ( 1 )  To  eliminate  B  as  a  neighbor  from  the  list(C). 

%  (2)  To  provide  a  pointer  to  A  in  the  list(C)  to  update  length(A,C)  and 
%  criterion(A,C)  in  list(C)  after  the  merging  of  A  and  B. 
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x=Start_N(C); 
while  x  ~=  0 
if  Reg_N(l,x)  ==  B 

Reg_N(l,x)=0;  %  Region  B  is  no  longer  a  neighbor  of  C. 

elseif  Reg_N(l,x)  =  A 

Pnt_to_A_inC=x;  %  Pointer  to  A  in  the  list(C). 

end 

x=Ptr_N(x);  %  Check  next  neighbor  to  see  if  it  is  B. 

end 


function  Reg_N  =  Criterion(A,C,pointer,M,Coeff,Reg_N) 

%  Already  Multichannel. 

%  To  update  the  criterion(A,C)  in  list(A). 

%  pointer  =  pointer  to  C  in  list(A). 

M_AC=M(A)+M(C); 

M_inv=l/M_AC; 

H=M(A)*M_inv*M(C); 

Reg_N(3,pointer)=Reg_N(2,pointer)-(Coeff(A)-Coeff(C))'*H*(Coeff(A)-Coeff(C)); 


function  Reg_N  =  Elim_and_Update(A,B>Coeff,M,Start_N,Ptr_N,Reg_N) 

%  After  taking  care  of  common  neighbors  and  before  merging  the  lists  of  A  - 
%  and  B  (B  will  be  absorbed  by  A.),  we  need  to  update  list(A)  and  list(B) 

%  as  follows. 

%  Updating  list(A): 

%  ( 1 )  Eliminate  B  as  neighbor  from  the  list(A). 

%  (2)  Update  criterion(A,C)  for  every  neighbor  C  in  list(A)  exept  neighb  B. 
%  (3)  Update  criterion(A,C)  in  list(C)  for  every  neighbor  C  in  list(A) 

%  exept  neighbor  B. 

%  Updating  list(B)  and  the  lists  of  the  members  of  list(B): 

%  (1 )  Eliminate  A  as  neighbor  from  the  list(B). 
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%  (2)  Update  criterion(A,C)  for  every  neighbor  C  in  list(B)  exept  neighb  A. 
%  (3)  Replace  B  by  A  in  the  list(C)  for  all  the  neighbors  C  of  B  in  list(B) 

%  except  neighbor  A  and  update  criterion(A,C)  in  list(C). 

%  NOTE  that  if  C  is  a  common  neighbor  of  A  and  B,  then  B  will  have  been 
%  replaced  by  0  at  the  common  neighbor  stage  (i.e.  function  "Elim_Pnt"). 
%  Thus  there  will  be  no  duplication  of  A  in  the  list(C). 

% . 


%  Updating  list(A). 

x=Start_N(A);  %  x=pointer  to  first  neighbor  C=Reg_N(l,x)  in  list(A). 

while  x  ~=  0  %  Ptr_N(x)=pnt  to  next  nghb  C=Reg_N(l,Ptr_N(x))  in  list(A). 

if  Reg_N(l,x)  =  B 

Reg_N(l,x)=0;  %Step(l). 

elseif  Reg_N(l,x)  ~=  0 

Reg_N=Criterion(A,Reg_N(  1  ,x),x,M,Coeff,Reg_N);  %Step  (2). 

Reg_N=Update_criterion(Reg_N(l  ,x),A,M,Coeff,Start_N,Ptr_N,Reg_N);  %Step  (3). 
end 

x=Ptr_N(x);  %  Check  next  neighbor  to  see  if  it  is  B. 

end 

% . 

%  Updating  list(B)  and  the  lists  of  the  members  of  list(B). 

x=Start_N(B);  %  x=pointer  to  first  neighbor  C=Reg_N(l ,x)  in  list(B). 

while  x  ~=  0 
if  Reg_N(l,x)  =  A 

Reg_N(  1  ,x)=0;  %  Step  ( 1 ). 

e  1  se  i  f  Reg_N(  1  ,x)  ~=  0 

Reg_N=Criterion(A,Reg_N(l,x),x,M,Coeff,Reg_N);  %  Step  (2). 

Reg_N=replaceB_by_A(A,B,Reg_N(  1  ,x),Coeff,M,Start_N,Ptr_N,Reg_N);  %  Step  (3). 
end 

x=Ptr_N(x);  %  Check  next  neighbor  to  see  if  it  is  A. 

end 


function  [Start_N,End_N,Ptr_N]  =  close_ranks(A,Start_N,End_N,Ptr_N,Reg_N) 
%  To  skip  over  the  addresses  in  the  pointer  Ptr_N  that  correspond  to 
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%  regions  in  list(A)  that  have  been  absorbed  and  adjust  Start_N  and  End_N. 

%  Note  that  when  region  "A"  absorbs  region  "B"  the  list  for  "A"  becomes  the 
%  union  of  list(A)  and  list(B),  and  it  has  at  least  two  new  zeros;  one  for 
%  "B"  in  list(A)  and  one  for  "A"  in  list(B).  Thus,  the  array  Reg_N(l,:)  has 
%  at  least  two  new  zeros  as  well. 

%  'The  lables  of  the  remaining  regions  are  the  non-zero  entries 
%  of  Reg_N(  1 

%  Here,  region  "A"  has  absorbed  region  "B";  however,  we  check  all  zeros  in 
%  list(A).  Extra  zeros  may  exist  in  list(A)  if  a  neighbor  of  "A"  has 
%  absorbed  a  region  that  is  also  a  neighbor  of  "A"  (see  "Elim_Pnt"). 

x=Start_N(A); 
if  x  ~=  0 
End_N(A)=x; 
while  Reg_N(l,x)  ==  0 
x=Ptr_N(x); 
if  x  ==  0 
return 
end 

Start_N(A):=x; 

End_N(A)=x; 
end 

while  Ptr_N(x)  ~=  0 
y=Ptr_N(x); 
if  y  ==  0 
return 
end 

while  Reg_N(l,y)  0 
Ptr_N(x)=Ptr_N(y); 
y=Ptr_N(y); 
if  y  ==  0 
return 
end 
end 
x=y; 

End_N(A)=y; 
end 
end 


%  Zeros  will  be  "skipped"  by  Ptr_N  in  list(A). 


%  Continue  looking  for  zeros  until  Ptr_N(x)=0. 
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function  [Start_B,End_B,Ptr_B,Reg_B]  =  New JBoundary(A,B, Start JB, End  JB,PtrJB,Reg_B) 


%  To  obtain  the  boundary  of  the  union  of  the  two  regions  A  and  B. 

%  The  addresses  of  the  boundary  points  of  A  are  Start_B(A),  Ptr_B(Start_B(A)), 

%  P tr__B( . . .Ptr_B(Start_B( A))),  ...End_B(A);  Ptr_B(End_B(A))=0  marks  the  end 
%  of  the  list.  Same  for  B. 

%  There  are  two  steps:  (1)  Delete  points  that  are  common  to  both  boundaries 
o/0  an(j  (2)  Concatenate  the  two  lists  and  close  ranks.. 

%  Step  1 . 

x=Start_B(A); 
y=Start_B(B); 
while  x  ~=  0 
while  y~=0 

if  Reg_B(x)  ==  Reg_B(y) 

Reg_B(x)=0; 

Reg  B(y)=0;  %  and/or  close  ranks! 

end 

y=Ptr_B(y); 

end 

x=Ptr_B(x); 

y=Start_B(B); 

end 

%  Step  2. 

Ptr_B(End_B(A))=Start_B(B);  %  To  concatenate  the  lists  for  A  and  B. 

End  B(A)=End_B(B);  %  The  end  of  the  list  for  A  is  now  at  the  end  of  B. 

Start  B(B)=0;  %  Region  B  has  been  absorbed  by  A. 

[Start  B,End_B,Ptr_B]=close_ranks(A,Start_B,End_B,Ptr_B,Reg_B); 


function  Reg_N  =  Update_criterion(C,A>M>Coeff,Start_N,Ptr_N,Reg_N) 

%  To  update  Criterion(C,A)  >n  list(C). 
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x=Start  N(C);  %  x  =  pointer  to  the  first  neighbor  of  C. 

while  x  ~=  0 

if  Reg_N(l,x)  =  A  %  Update  Reg_N(3,x)=criterion(C,A). 

Reg_N=Criterion(C,A>x,M,Coeff,Reg_N); 

end 

x=Ptr  N(x);  %  Check  next  neighbor  to  see  if  it  is  A. 

end 


function  Reg_N  =  replaceBJby_A(A,B,C,Coeff,M,Start_N,Ptr_N,Reg_N) 


%  To  replase  B  by  A  in  list(C)  and  update  criterion(A,C)  in  list(C). 

x=Start_N(C);  %  x  =  pointer  to  the  first  neighbor  of  C. 

while  x  ~=  0 

if  Reg_N(  1  ,x)  ==  B  %  Replace  B  by  A  and  update  criterion(A,C)  in 

%  list(C). 

Reg_N(  1  ,x)=A;  %  Region  B  has  been  replaced  by  A. 

Reg_N=Criterion(C,A,x,M,Coeff,Reg_N); 
end 

x=Ptr__N(x);  %  Check  next  neighbor  to  see  if  it  is  B. 

end 


function  [Start, End, Pointer]  =  merge_pix(A,B,Start,End,Pointer) 

%  To  update  the  pointer  "Pointer"  when  a  merging  of  regions  A  and  B  occurs. 
%  The  pointer  points  to  the  pixels  in  regions  A  and  B.  The  new  region  is 
%  labled  A. 

Pointer(End(A))=Start(B); 

End(A)=End(B); 

Start(B)=0; 


function  [Strt_S,End_S,P_S]  =  cls_rnks(Strt_S,End_S,P_S, Start) 

%  To  skip  over  the  addresses  in  the  pointer  P_S  that  correspond  to 
%  a  region  that  has  been  absorbed  and  adjust  Strt_S  and  End_S. 
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%  Note  that  when  region  "A"  absorbs  region  "B"  the  pixels  for  "A"  becomes 
%  the  union  of  pixels(A)  and  pixels(B),  and  Start  has  a  new  zero  at 
%  Start(B). 


x~St:rt_S; 

End_S=x; 
while  Start(x)  ==  0 
x=P_S(x); 
if  x  ==  0 
return 
end 

Strt_S=x; 

End_S=x; 

end 

while  P_S(x)  ~=  0 
y=P_S(x); 

if  y  ==  0,  return,  end 
while  Start(y)  =  0 

P_S(x)=P_S(y);  %  A  zero  will  be  "skipped"  by  P_S. 

y=P_S(y); 
if  y  ==  0 
return 
end 
end 

x=y;  %  Check  next  entry. 

End_S=y; 

end 


PART  IV 

COLLECTING  REGIONS  &  BOUNDARIES  &  THE  PC-SEGMENTED  IMAGE 

function  [Regns, Segmentation, FinaI_Regions]  =  Display_Reg_Pix(Start_S,Ptr_S,Start,... 

Reg_Pix,C) 

%  To  collect  the  lables  of  all  the  regions  in  the  final  segmentation  and 
%  put  them  in  the  array  Final_Regions.  The  size  of  Final_Regions  is  the 
%  number  of  regions  in  the  segmentation. 
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%  To  assign  the  lable  "A"  to  every  pixel  in  region  A  for  every  region  A  in 
%  the  final  segmentation  and  put  it  in  the  lxnm-array  Regns. 

%  To  assign  the  piecewise  constant  value  C(A)  to  every  pixel  in  region  A 
%  for  every  region  A  in  the  final  segmentation  and  put  it  in  the  lxnm-array 
%  Segmentation. 

%  Start_S  and  Ptr_S  are  associated  with  the  lables  of  regions  . 

%  Start  and  RegPix  are  associated  with  the  pixels  ("start  &  pointer"). 


Final_Regions=[]; 


A=Start_S; 
while  A  ~=  0 
if  Start(A)  ~=  0 

Final_Regions=[Final_Regions  A]; 
x=Start(A); 

Regns(x)=A; 

Segmentation(x)=C(A); 

%  Seginentation(:,x)=C(:,A); 
x=Reg_Pix(x); 
while  x  ~=  0 
Regns(x)=A; 
Segmentation(x)=C(A); 
x=Reg_Pix(x); 
end 
end 

A=Ptr_S(A); 


%  First  region  in  the  final  segmetation  has  lable  A. 


%  A  is  a  final  region. 

%  First  pixel  in  region  A. 

%  Pixel  x  receives  lable  A. 

%  In  the  segmented  image,  pixel  x  has  PC-value  C(A). 
%  MULTICHANNEL  PC-value. 

%  Next  pixel  in  region  A. 

%  All  pixels  x  in  A  receive  lable  A. 

%  All  pixels  x  in  A  have  PC-value  C(A). 

%  Next  pixel  in  region  A. 


%  Next  region  in  the  final  segmentation. 


function  Image  =  Vect_to_Img_l(m,V) 

%  This  function  transforms  a  vector  representation  of  an  image  V  into  a 
%  matrix  representation  "Image".  The  matrix  "Image"  is  nxm. 

%  Jorge  M.  Martin,  NAWC-WPNS  Code  474400D,  China  Lake,  CA.  April  1996. 

L=length(V);  %  L=nm  should  be  a  multiple  of  m. 

for  k=l  :L 
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I  mage(fix((k- 1  )/m)+ 1  ,rem(k- 1  ,m)+ 1  )=V  (k); 
end 


function  ith_region  =  Display_ith_reg(A,Start,Reg_Pix,nm,C) 

%  A  is  the  lable  of  the  ith-region.  The  nm-row  vector  "ith_region"  will 
%  contain  zeros  everywhere  exept  at  locations  corresponding  to  the  pixels 
%  of  region  A.  At  these  locations  it  will  either  contain  the  lable  "A"  or 
%  the  constant  value  C(A)  of  the  approximation  to  the  original  image  that 
%  is  associated  with  region  A  of  the  segmented  image. 


ith_region=zeros(  1  ,nm); 


if  Start(A)  ~=  0 
x=Start(A); 
ith_region(x)=A; 

%  ith_region(x)=C(A); 
x=Reg_Pix(x); 
while  x  ~=  0 

ith_region(x)=A; 

%  ith_region(x)=C(A); 

x=Reg_Pix(x); 

end 

end 


%  x  is  the  First  pixel  in  region  A. 

%  Lable  A  put  at  location  x...,or... 

%  constant  value  C(A)  put  at  location  x. 
%  Next  pixel  in  region  A. 

%  Lable  A  put  at  location  x...,or... 

%  constant  value  C(A)  put  at  location  x. 
%  Next  pixel  in  region  A. 


function  [Bndrs,Reg_B]  =  Display_Boundaries(Final_Regions,Start_B,Ptr_B,Reg_B,nm) 

%  To  put  a  one  on  every  k  in  { l,2,...,nm}  that  is  a  boundary  point. 

Bndrs=zeros(l,nm); 

N=length(Final_Regions); 
for  i=l  :N 

x=Start_B(Final_Regions(i)); 
while  x  ~=  0 
if  Reg_B(x)  >  nm 
Reg_B(x)=Reg_B(x)-nm; 
end 

if  Reg_B(x)  ~=  0 
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Bndrs(Reg_B(x))=l ; 
end 

x=Ptr_B(x); 

end 


function  ith_Boundary  =  DisplJthJBndry(A, Start JB,Ptr_B,Reg_B,nm) 

%  To  get  the  pixels  of  the  ith-boundary. 

ith_Boundary=zeros(  1  ,nm); 
x=Start_B(A); 
while  x  ~=  0 
if  Reg_B(x)  ~=  0 
ith_Boundary(Reg_B(x))=  1 ; 
end 

x=Ptr_B(x); 

end 
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